module hcncool_module implicit none ! ------------------------------------------------------------- ! Attempt at incorporating the hcncooling code into the ! Titan TGCM Framework. ! Jared Bell 02-23-05 ! ------------------------------------------------------------- INTEGER, PARAMETER :: NLINES = 65 REAL,DIMENSION(NLINES) :: VI REAL,DIMENSION(NLINES) :: VJ REAL,DIMENSION(NLINES) :: GAMMA REAL,DIMENSION(NLINES,9) :: P ! PRAT is the Cooling Ratio Coefficients from Matlab 03-29-05 REAL,DIMENSION(27,4) :: PRAT * ----------------------------------------------------------------- * FUNDAMENTAL CONSTANTS-------------------------------------------- * ----------------------------------------------------------------- * C = Speed of light in m/s * C2 = Speed of light in cm/s (useful for cm^-1 to Hz conversions) * KB = Boltzmann's Constant * H = Planck Constant * AMU = Atomic Mass Unit in kg * MHCN = Mass of HCN Molecule in kg * RT = Titan Radius in km * MT = Titan Mass in kg * G = Fundamental Gravitational Constant in SI units * ----------------------------------------------------------------- * ----------------------------------------------------------------- REAL, PARAMETER :: C = 2.99792D8 REAL, PARAMETER :: C2 = 2.99792D10 REAL, PARAMETER :: KB = 1.38065D-23 REAL, PARAMETER :: H = 6.62607D-34 REAL, PARAMETER :: PI = 3.141592653589793D0 REAL, PARAMETER :: AMU = 1.6605D-27 REAL, PARAMETER :: MHCN = 27.0D0*AMU REAL, PARAMETER :: RT = 2.575D3 REAL, PARAMETER :: MT = 1.3500336D3 REAL, PARAMETER :: G = 6.6726D-11 REAL, PARAMETER :: ATMTOPAS = 1.0D0/(1.01295D5) contains * --------------------------------------------------------------------- SUBROUTINE hcncool(TI,NCH4I,NH2I,NHCNI,barmI,fcpI,ZI, | xnmbari,TCOOL,blev,elev,blon,elon,lat) * --------------------------------------------------------------------- * Purpose: * ======= * This is the Cooling Routine to be used for the Titan TGCM * This version incorporates all heating and cooling terms into * the full line by line radiative transfer calculations. * Date Programmer Version * ========== ================== ============== * 09-08-04 Jared Bell Original(1.0) * 12-02-04 Jared Bell Update (1.1) * 12-07-04 Jared Bell Update (1.2) * 01-27-05 Jared Bell Consistency Check * 01-27-05 Jared Bell Added Heating 1 * 01-27-05 Jared Bell Added Heating 2 * 01-28-05 Jared Bell Added Dynamic Memory * ----------------------------------------------------------------- * INPUTS: * ====== * Essentially, this code requires 8 inputs from the main * code: Temperature, HCN Density, and Total Density, Altitude * ----------------------------------------------------------------- * Currently, the code simply extracts these three parameters from * a static background profile from the Yung, Allen, and Pinto paper * on Titan (1984). * ------------------------------------------------------------------ * BACKGROUND ON CODE: * ================== * Essentially, the code assumes the INPUTS in a paticular format. * Since there are 65 Rotational Lines (Denoted NLINES in the CODE) * and 27 Pressure levels (Denoted NLVLS in the code -- * The code expects that the input be given on cubes of size * NLVLS x NLON * ------------------------------------------------------------------ * INPUT FORMAT: * ============ * Currently, the code makes ample use of the Element by Element Matrix * Mathematics enabled in FORTRAN 90. Thus all Inputs are currently assumed * to be matrices of size (NLVLS,NLON). * ------------------------------------------------------------------ * DECLARATIONS SECTION --------------------------------------------- * ------------------------------------------------------------------ * NLINES, NLVLS, are the number of Rotational Lines and Height/Pressure * Levels Respectively. * P is an array carrying Polynomial Coefficients for the Intensity * Interpolation. * ----------------------------------------------------------------- * * ----------------------------------------------------------------- * VARIABLE DECLARATIONS-------------------------------------------- * ----------------------------------------------------------------- use cons_module, only: rmassinv_ch4, rmassinv_hcn, | rmassinv_h2 IMPLICIT NONE * INPUTS:---------------------------------------------------------- ! ----------------------------------------------------------------- ! These are 2-D arrays that are passed as inputs from dynamics.F ! They are initialized with the sizes passed from dynamics.F with ! Dimension(lev0:lev1,lon0:lon1). ! INTEGER,INTENT(IN):: blon,elon,blev,elev,lat REAL,INTENT(IN),DIMENSION(blev:elev,blon-2:elon+2) :: | TI, | NCH4I, | NH2I, | NHCNI, | barmI, | fcpI, | ZI, | xnmbari ! * ----------------------------------------------------------------- * OUTPUTS:--------------------------------------------------------- * REAL,INTENT(OUT),DIMENSION(blev:elev,blon-2:elon+2) :: | TCOOL * ----------------------------------------------------------------- * LOCAL VARS:------------------------------------------------------ INTEGER :: i,j,l,m,k,n,ii INTEGER :: NLVLS,NLON INTEGER :: | start, | end, | rate INTEGER,DIMENSION(2) :: | index logical, parameter :: | debug = .true. logical :: | done = .false. REAL,DIMENSION(blev:elev,blon-2:elon+2) :: | TCOOLa REAL :: | FTH, | TIME REAL,DIMENSION(NLINES,blev:elev) :: | S, | INTENSITY REAL,DIMENSION(elon-blon+1,blev:elev) :: | T, | NCH4, | NH2, | NHCN, | barm, | xnmbar, | rho, | fcp, | Z, | TCOOL2, | TCOOL2a, | COOL, | COOL2, | HEAT1, | HEAT2, | SMCOOL, | SMHEAT1, | SMHEAT2, | RATIO REAL,DIMENSION(elon-blon+1,NLINES,blev:elev) :: | PRODCOOL, | PRODHEAT1, | PRODHEAT2, | gcool, | gheat1, | gheat2, | INTEN, | ALPHAD, | DELTAV, | VULIM, | VLLIM, | VJM, | TM, | ZM, | VTH, | NHCNM REAL,DIMENSION(blev:elev,blon-2:elon+2) :: | NHCNT NLVLS = elev - blev + 1 NLON = elon - blon + 1 HEAT1 = 0.0D0 HEAT2 = 0.0D0 ! Transfer inputs to local arrays: do i = blon,elon ii = i - blon + 1 T(ii,:) = TI(:,i) NCH4(ii,:) = NCH4I(:,i)*xnmbari(:,i)*rmassinv_ch4 NH2(ii,:) = NH2I(:,i)*xnmbari(:,i)*rmassinv_h2 NHCN(ii,:) = NHCNI(:,i)*xnmbari(:,i)*rmassinv_hcn ! Number density (cm^-3) barm(ii,:) = barmI(:,i) xnmbar(ii,:) = xnmbari(:,i) fcp(ii,:) = fcpI(:,i) Z(ii,:) = ZI(:,i)/1.0D4 end do rho = (xnmbar*AMU*1.0D3) * ----------------------------------------------------------------- * END ALL DECLARATIONS -------------------------------------------- * ----------------------------------------------------------------- * * ------------------------------------------------------------------------ * INTERPOLATING INTENSITY DATA TO TEMPERATURE GRID ----------------------- * ------------------------------------------------------------------------ ! Define s(nlines,blev:elev), and inten(nlon,nlines,nlev) DO j = 1, NLINES DO i = 1, NLON DO k = 1, NLVLS S(j,k) = P(j,9) + P(j,8)*T(i,k) + P(j,7)*(T(i,k)**2.0D0) | + P(j,6)*(T(i,k)**3.0D0)+P(j,5)*(T(i,k)**4.0D0) | + P(j,4)*(T(i,k)**5.0D0)+P(j,3)*(T(i,k)**6.0D0) | + P(j,2)*(T(i,k)**7.0D0)+P(j,1)*(T(i,k)**8.0D0) INTEN(i,j,k) = C2*(1.0D-4)*S(j,k) END DO END DO END DO * ------------------------------------------------------------------------ * END INTERPOLATING INTENSITY DATA TO TEMPERATURE GRID ------------------- * ------------------------------------------------------------------------ * DO k = 1, NLVLS RATIO(:,k) = PRAT(k,4) + | PRAT(k,3)*T(:,k) + | PRAT(k,2)*(T(:,k)**2.0D0) + | PRAT(k,1)*(T(:,k)**3.0D0) END DO * * ------------------------------------------------------------------------ * BUILDING FULL 3-D MATRICES TM,VJM,ZM,NHCNM * ------------------------------------------------------------------------ * --------------- * FOR LOOP BEGIN- * --------------- DO k = 1, NLVLS DO i = 1, NLON VJM(i,:,k) = VJ(:) ZM(i,:,k) = Z(i,k) TM(i,:,k) = T(i,k) NHCNM(i,:,k) = NHCN(i,k) END DO END DO * --------------- * END FOR LOOP -- * --------------- * ------------------------------------------------------------------------ * END BUILDING 3-D MATRICES ZM,TM,VJM,NHCNM * ------------------------------------------------------------------------ * * * ------------------------------------------------------------------------ * DOPPLER HALFWIDTHS AT GIVEN TEMPERATURES ------------------------------- * ------------------------------------------------------------------------ * FTH = SQRT(2.0D0*(KB/MHCN)) VTH = FTH*SQRT(TM) ALPHAD = (VJM*(VTH/C)) DELTAV = 2.75D0*ALPHAD VULIM = VJM + DELTAV VLLIM = VJM - DELTAV * ------------------------------------------------------------------------ * END DOPPLER HALFWIDTHS AT GIVEN TEMPERATURES --------------------------- * ------------------------------------------------------------------------ * * ------------------------------------------------------------------------ * Cooling CALCULATIONS * ------------------------------------------------------------------------ ! TIME = 0.0D0 ! CALL system_clock(COUNT = start, COUNT_RATE = rate) CALL fgausscool(VLLIM,VULIM,ALPHAD,TM,VJM,gcool,NLVLS,NLON) ! PRODCOOL = gcool*INTEN*(4.0D0*PI) PRODCOOL = gcool*INTEN*(2.0D0*PI) SMCOOL = SUM(PRODCOOL,2) ! COOL = SMCOOL*NHCN*(1.0D7)*(1.0D0/(rho*fcp))*(3600.*24.) !(K/Day) COOL = SMCOOL*NHCN*(1.0D7)*(1.0D0/(rho*fcp)) ! K/sec COOL2 = SMCOOL*NHCN*(1.0D7) ! ergs/(cm^3 s) ! The factor of 1/2 scales the total cooling to cool-to-space values TCOOL2 = RATIO*COOL ! This scales cooling to K/s or K/day TCOOL2a = RATIO*COOL2 ! This scales cooling to erg/cm3/sec * ------------------------------------------------------------------------ * Cooling CALCULATIONS * ------------------------------------------------------------------------ !--------------------------------------------------------------------------- DO i = blon,elon ii = i - blon + 1 DO k = 1, NLVLS TCOOL(k,i) = TCOOL2(ii,k) TCOOLa(k,i) = TCOOL2a(ii,k) NHCNT(k,i) = NHCN(ii,k) END DO END DO ! ! call addfsech('TCOOL', ' ' , ' ' ,TCOOL(:,blon:elon),blon, ! ! elon,NLVLS,NLVLS,lat) ! call addfsech('TCOOLa', ' ' , ' ' ,TCOOLa(:,blon:elon),blon, ! ! elon,NLVLS,NLVLS,lat) END SUBROUTINE hcncool *--------------------------------------------------------------------------- * SUBROUTINE DECLARATION SECTION *--------------------------------------------------------------------------- *--------------------------------------------------------------------------- * FGAUSSCOOL DECLARATION *--------------------------------------------------------------------------- ! SUBROUTINE fgausscool(A,B,ALPHAD,TM,VJM,gcool,NLVLS,NLON) IMPLICIT NONE ! INPUT SECTION --------------------------------------------------------- INTEGER, INTENT(IN) :: NLVLS,NLON REAL,INTENT(IN), DIMENSION(1:NLON,1:NLINES,1:NLVLS) :: | A, | B, | ALPHAD, | TM, | VJM ! REAL,INTENT(IN),DIMENSION(NLON,NLINES,NLVLS):: ALPHAL ! Outputs ----------------------------------------------------- REAL,INTENT(OUT),DIMENSION(1:NLON,1:NLINES,1:NLVLS):: gcool ! LOCAL VARIABLES SECTION ----------------------------------------------- INTEGER :: i,j,k,m,n,status REAL,DIMENSION(1:NLON,1:NLINES,1:NLVLS) :: | S1, | OLD, | X1, | X2, | Y1, | Y2, | Y REAL,DIMENSION(8,4) :: D REAL,DIMENSION(8,4) :: F ! ! END LOCAL VARIABLES SECTION -------------------------------------------- n = 16 S1 = 0.0D0 D(1,1) = 1.0D0 D(1:2,2) = (/.6521451548D0, .3478548451D0/) D(1:4,3) = (/.3626837833D0, .3137066458D0, .2223810344D0, ! .1012285362D0/) D(1:8,4) = (/.1894506104D0, .1826034150D0, .1691565193D0, ! .1495959888D0, ! .1246289712D0, .0951585116D0, .0622535239D0, .0271524594D0/) F(1,1) = .5773502691D0 F(1:2,2) = (/.3399810435D0, .8611363115D0/) F(1:4,3) = (/.1834346424D0, .5255324099D0, .7966664774D0, ! .9602898564D0/) F(1:8,4) = (/.0950125098D0, .2816035507D0, .4580167776D0, ! .6178762444D0, & .7554044084D0, .8656312023D0, .9445750230D0, .9894009350D0/) j = 1 *---------------- * WHILE LOOP BEGIN *---------------- DO IF (2**j .EQ. n) EXIT j = j + 1 END DO *---------------- * END WHILE LOOP *--------------- * * --------------- * FOR LOOP BEGIN- * --------------- DO k = 1 , n/2 X1 = (F(k,j)*(B - A) + A + B)/2.0D0 X2 = (-F(k,j)*(B - A) + A + B)/2.0D0 CALL coolgauss(X1,ALPHAD,TM,VJM,Y1,NLVLS,NLON) CALL coolgauss(X2,ALPHAD,TM,VJM,Y2,NLVLS,NLON) Y = Y1 + Y2 OLD = S1 S1 = OLD + D(k,j)*Y END DO * --------------- * END FOR LOOP -- * --------------- gcool = (B-A)*(S1/2) END SUBROUTINE fgausscool * *--------------------------------------------------------------------------- * COOLGAUSS DECLARATION *--------------------------------------------------------------------------- SUBROUTINE coolgauss(V,ALPHAD,TM,VJM,ans,NLVLS,NLON) IMPLICIT NONE * * INPUT SECTION ------------------------------------------------------------ INTEGER, INTENT(IN) :: NLVLS,NLON REAL,INTENT(IN),DIMENSION(NLON,NLINES,NLVLS)::V REAL,INTENT(IN),DIMENSION(NLON,NLINES,NLVLS)::ALPHAD REAL,INTENT(IN),DIMENSION(NLON,NLINES,NLVLS)::TM REAL,INTENT(IN),DIMENSION(NLON,NLINES,NLVLS)::VJM * REAL,INTENT(IN),DIMENSION(NLON,NLINES,NLVLS)::ALPHAL REAL,INTENT(OUT),DIMENSION(NLON,NLINES,NLVLS)::ans * * LOCAL VARIABLES --------------------------------------------------------- INTEGER status,i,j,m REAL,DIMENSION(NLON,NLINES,NLVLS) :: PHI REAL,DIMENSION(NLON,NLINES,NLVLS) :: ARG REAL,DIMENSION(NLON,NLINES,NLVLS) :: ALPHA REAL,DIMENSION(NLON,NLINES,NLVLS) :: B REAL,DIMENSION(NLON,NLINES,NLVLS) :: FACTOR REAL,DIMENSION(NLON,NLINES,NLVLS) :: TOTAL * VARIABLE SECTION --------------------------------------------------------- * ALPHA = (H/KB)*(V/TM) B = ((2.0D0*H)/(C**2))*(V**3)*(1.0D0/(EXP(ALPHA) - 1.0D0)) FACTOR = 1.0D0/(SQRT(PI)*ALPHAD) ARG = ((V - VJM)/ALPHAD)**2 PHI = FACTOR*EXP(-ARG) TOTAL = PHI*B ans = TOTAL END SUBROUTINE coolgauss ! !-------------------------------------------------------------------- ! subroutine init_hcn ! This is for titan tgcm radiation code (see hao:/home/tgcm/titan/rad) implicit none INTEGER :: i,j,m,l integer, parameter :: np = 9 integer, parameter :: nlevs = 27 integer, parameter :: nr = 4 integer :: lu=2 call read_nc(VI,GAMMA,P,PRAT,NLINES,nlevs,np,nr) DO j = 1, NLINES VJ(j) = C2*VI(j) END DO end subroutine init_hcn !----------------------------------------------------------------------- subroutine read_nc(vi,gamma,p,prat,nlines,nlevs,np,nr) ! ! Example routine to read netcdf file ncfile (in module data above), ! written by sub write_nc. This code can be used in the titan hcncooling ! module. ! implicit none #include ! ! Args: integer, intent(in) :: nlines,np,nr,nlevs real, intent(in), dimension(nlines) :: vi real, intent(in), dimension(nlines) :: gamma real, intent(in), dimension(nlines,np) :: p real, intent(in), dimension(nlevs,nr) :: prat ! Local: character(len=80) :: ncfile = 'hcncool.nc' integer :: i,ncid,istat integer :: nlevrd,id_nlines,nlinesrd,id_np,id_nlevs,id_nr, | nprd,id_vi,id_gamma,id_p,id_prat character(len=80) :: description,create_date character(len=120) :: msg ! ! Open existing netcdf file: istat = nf_open(ncfile,NF_NOWRITE,ncid) if (istat /= NF_NOERR) then write(msg,"('Error opening netcdf file ',a)") trim(ncfile) call handle_ncerr(istat,trim(msg)) endif write(6,"(/,'Opened netcdf file ',a)") trim(ncfile) ! ! Get global attributes: description = ' ' istat = nf_get_att_text(ncid,NF_GLOBAL,"Description",description) ! write(6,"('Description: ',a)") trim(description) ! create_date = ' ' istat = nf_get_att_text(ncid,NF_GLOBAL,"Creation_date", | create_date) ! write(6,"('Create_date: ',a)") trim(create_date) ! ! Read rotational line parameters: ! ! VI: istat = nf_inq_varid(ncid,'VI',id_vi) if (istat /= NF_NOERR) | call handle_ncerr(istat,"Error getting VI var id") istat = nf_get_var_double(ncid,id_vi,vi) if (istat /= NF_NOERR) | call handle_ncerr(istat,"Error getting VI var") ! write(6,"('VI=',/,(6e12.4))") vi ! ! GAMMA: istat = nf_inq_varid(ncid,'GAMMA',id_gamma) if (istat /= NF_NOERR) | call handle_ncerr(istat,"Error getting GAMMA var id") istat = nf_get_var_double(ncid,id_gamma,gamma) if (istat /= NF_NOERR) | call handle_ncerr(istat,"Error getting GAMMA var") ! write(6,"('GAMMA=',/,(6e12.4))") gamma ! ! Read intensity coefficients: ! ! P: istat = nf_inq_varid(ncid,'P',id_p) if (istat /= NF_NOERR) | call handle_ncerr(istat,"Error getting P var id") istat = nf_get_var_double(ncid,id_p,p) if (istat /= NF_NOERR) | call handle_ncerr(istat,"Error getting P var") ! do i=1,np ! write(6,"('i=',i3,' P(:,i)=',/,(6e12.4))") i,p(:,i) ! enddo ! ! Read RATIO coefficients: ! ! PRAT: istat = nf_inq_varid(ncid,'PRAT',id_prat) if (istat /= NF_NOERR) | call handle_ncerr(istat,"Error getting PRAT var id") istat = nf_get_var_double(ncid,id_prat,prat) if (istat /= NF_NOERR) | call handle_ncerr(istat,"Error getting PRAT var") ! do i=1,nr ! write(6,"('i=',i3,' PRAT(:,i)=',/,(6e12.4))") i,prat(:,i) ! enddo ! end subroutine read_nc !----------------------------------------------------------------------- subroutine handle_ncerr(istat,msg) ! ! Handle a netcdf lib error: ! integer,intent(in) :: istat character(len=*),intent(in) :: msg ! write(6,"(/72('-'))") write(6,"('>>> Error from netcdf library:')") write(6,"(a)") trim(msg) write(6,"('istat=',i5)") istat write(6,"(72('-')/)") ! return stop end subroutine handle_ncerr !----------------------------------------------------------------------- !-----END OF HCN COOLING--------------------------------------- !----------------------------------------------------------------------- end module hcncool_module