#include "dims.h" ! ! am_02/02: DT1DTS added (DT0DTS / abs[sin(Im)]) ! DT0DTS calculated from DT1DTS (D(THETA0)/D(THETAS)) ! module cons_module implicit none ! ! Define model constants. ! Parameter constants are cons_module module data and are accessed ! in subroutines via use-association. ! Derived constants are cons_module module data, and are calculated ! in sub init_cons (contained in cons_module). ! Sub init_cons is called by sub init (init_mod.F). ! Parameters referenced here are in params.h (s.a., dims.h) ! #include "params.h" ! ! Integer parameter constants: ! integer,parameter :: | imax=zimx, ! number of longitudes | jmax=zjmx, ! number of latitudes | kmax=zkmx, ! number of levels (midpoints) | imaxh=imax/2, | imaxp2=imax+2, | imaxp4=imax+4, | jmaxm1=jmax-1, | jmaxm2=jmax-2, | jmaxh=jmax/2, | js=-jmaxh, | kmaxm1=kmax-1, | kmaxp1=kmax+1, | last=imaxp4*kmax-2, | len1=imaxp4, | len2=imaxp4*kmax, | len3=imaxp4*kmaxp1, | long=len2-4 ! ! Real parameter constants: ! real,parameter :: | re = 6.37122e8, ! earth radius C(51) | re_inv = 1./re, ! inverse of earth radius C(52) | avo = 6.023e23, ! avogadro number C(85) | boltz = 1.38E-16, ! boltzman's constant C(84) | p0 = 5.0e-4, ! standard pressure C(81) | gask = 8.314e7, ! gas constant C(57) | dipmin = .005, ! minimum mag dip angle (tiegcm) ! | dipmin = 0.17, ! minimum mag dip angle (timegcm) | brn2d = 0.6, ! ! | tbound = 233.7, ! background tn at lower boundary (timegcm) | tbound = 181.0, ! background tn at lower boundary (tiegcm) | atm_amu = 28.9, ! mean mass of surface atmosphere C(24) | shapiro = 3.0E-2, ! shapiro smoother constant C(26) | dtsmooth = 0.95, ! time smoothing constant C(30) | dtsmooth_div2 = 0.5*(1.-dtsmooth), ! C(31) ! | tgrad = 20., ! TN gradient from old /RUNMDS/ (timegcm) | tgrad = 6., ! TN gradient from old /RUNMDS/ (tiegcm) | nob(zjmx) = 4.e6, ! N(NO) LBC from old /RUNMDS/ | avto = 4.0e-12, ! from old /RUNMDS/ (not used) ! | hor(zjmx) = 1., ! horizontal eddy diffusion (timegcm) | hor(zjmx) = .25, ! horizontal variation of eddy diffusion and ! eddy thrmal conductivity. If unity, value ! of KE at bottom is 5.e-6 (old /RUNMDS/) | prndtl = 1., ! prandtl number | evergs = 1.602e-12, ! 1 eV = 1.602e-12 ergs | tsurplus=5.11*evergs ! surplus heat per event (ergs) C(45) ! ! Many expressions require x/rmass, but its more efficient on many ! platforms to multiply rather than divide, so set rmassinv = 1./rmass ! here, and use x*rmassinv in the code. ! real,parameter :: | rmass(3) = (/32.,16.,28./), ! o2,o,n2 | rmass_o2 = 32., rmass_o = 16., rmass_n2 = 28., rmass_o3 = 48., | rmassinv_o2 = 1./rmass_o2, | rmassinv_o = 1./rmass_o, | rmassinv_n2 = 1./rmass_n2, | rmassinv_o3 = 1./rmass_o3 ! ! 2/00: these were in modsrc.snoe (tgcm13mt), but were unused. ! Low-energy protons: ! real,parameter :: ! alfalp = 10., ! efluxlp = 1.e-20 ! ! Model derived constants (see sub init_cons in this module): ! real :: | pi, ! set with 4*atan(1) C(110) | dphi, ! delta lat (pi/nlat) C(2) | dphi_2div3, ! 2./(3.*dphi) C(12) | dphi_1div12, ! 1./(12.*dphi) C(13) | dlamda, ! delta lon (2pi/nlon) C(1) | dlamda_2div3, ! 2./(3.*dlamda) C(10) | dlamda_1div12, ! 1./(12.*dlamda) C(11) | dt, ! time step C(4) | dtx2, ! 2*dt C(6) | dtx2inv, ! 1./(2*dt) C(7) | freq_3m3, ! frequency of 2-day wave (rad/sec) C(21) | freq_semidi, ! frequency of semidiurnal tide (rad/sec) C(23) | freq_ann, ! frequency of annual tide C(25) | expz(zkmxp), ! exp(-z) at midpoints | expzmid, ! exp(-.5*dz) C(86) | expzmid_inv, ! 1./expzmid C(87) | rmassinv(3), ! inverse of rmass | t0(zkmxp), ! set by sub lowbound (bndry_mod.F) | cs(-1:zjmx+2), ! cos(phi) | sn(zjmx), ! sin(phi) | tn(zjmx), ! tan(phi) | cor(zjmx), | grav, ! accel due to gravity (dependent on lower boundary) | dzgrav, ! grav/gask C(65) | difk(zkmxp), ! background eddy diffusion | dift(zkmxp), ! background thermal conductivity | xmue(zkmxp) ! eddy viscosity (?) ! ! kut for timegcm: ! integer,parameter :: kut(zjmx) = ! | (/1,2,3,5,6,8,11,16,21,26,31,36,36,36,36,36,36,36,36,36,36, ! | 36,36,36,36,31,26,21,16,11,8,6,5,3,2,1/) ! ! kut for tiegcm: integer,parameter :: kut(zjmx) = | (/1,2,3,5,6,7,9,10,11,13,14,15,17,17,17,17,17,17,17,17,17,17,17, | 17,15,14,13,11,10,9,7,6,5,3,2,1/) ! integer,parameter :: kflds(nflds) = (/ | 3 ,3 ,3 ,3 ,3 ,3 , | 3 ,3 ,3 ,3 ,3 ,3 , | 3 ,3 ,3 ,3 ,1 ,3 , | 3 ,3 ,3 ,3 ,3 ,3 , | 3 ,3 ,3 ,2 ,2 ,2 , | 2 ,2/) ! integer,parameter :: kpole(nflds) = (/ | 1 ,-1 ,-1 ,1 ,1 ,1 , | 1 ,1 ,0 ,1 ,1 ,0 , | 0 ,1 ,1 ,1 ,1 ,1 , |-1 ,-1 ,1 ,1 ,1 ,1 , | 1 ,1 ,1 ,0 ,0 ,0 , | 0 ,0/) ! character(len=8) :: nflds_lab(nflds) = | (/'NT ','NU ','NV ','NPS ','NPS2 ', | 'NPN4S ','NPNO ','NOP ','NPN2D ','NTI ', | 'NTE ','NE ','NO2P ','NW ','NZ ', | 'NPHI ','NPHIH ', | 'NTNM ','NUNM ','NVNM ','NPSNM ','NPS2NM ', | 'NN4SNM ','NPNONM ','NOPNM ','NMS ','NVC ', | 'NKLDU ','NKLDV ','NKLDT ','NKLDPS ','NKLDP2 '/) ! character(len=8) :: nphys_lab(nphys) = | (/'NRH ','NFLH ','NFPH ','NQDH ','NKMH ', | 'NPSDH ','NPSDH2 ','NLXX ','NLYY ','NLXY ', | 'NLYX ','NKT ','NKM ','NCP ','NWTI ', | 'NWTE ','NNO2 ','NNO ','NNN2 ','NPDUM ', | 'NNVO2 ','NNVO ','NNVN2 ','NQO2P ','NQOP ', | 'NQN2P ','NQNOP ','NQNP ','NQOP2D ','NQOP2P ', | 'NQTEF ','NN2P ','NNP ','NNOP ','NQ ', | 'NRJ ','NUI ','NVI ','NWI '/) ! integer :: ncols ! set by init_buff ! contains !----------------------------------------------------------------------- subroutine init_cons use input_module,only: step ! ! Set derived constants (this is called from sub init in init_module) ! ! Local: real :: z,expdz,phi real :: omega = 7.292E-5 integer :: i,k,j,nfsech,js ! pi = 4.*atan(1.) ! C(110) dphi = pi/float(zjmx) ! C(2) dphi_2div3 = 2./(3.*dphi) ! C(12) dphi_1div12 = 1./(12.*dphi) ! C(13) dlamda = 2.*pi/float(zimx) ! C(1) dlamda_2div3 = 2./(3.*dlamda) ! C(10) dlamda_1div12 = 1./(12.*dlamda) ! C(11) ! ! expz(kmax) is exp(-zp) at midpoints: ! expz (will replace EXPS) (expz(zkmxp) not used). ! expz(:) = 0. ! init ! ! bottom midpoint z = zsb + 1/2 deltaz (deltaz==dz==0.5 or 0.25) ! (zsb and dz are in params.h) z = zsb+.5*dz expz(1) = exp(-z) expdz = exp(-dz) difk(1) = 5.0e-6 dift(1) = 5.0e-6/prndtl xmue(1) = 5.0e-6 do k=2,zkmx expz(k) = expz(k-1)*expdz difk(k) = difk(k-1)*expdz xmue(k) = difk(k) dift(k) = dift(k-1)*expdz enddo difk(zkmxp) = difk(zkmx)*expdz dift(zkmxp) = dift(zkmx)*expdz xmue(zkmxp) = difk(zkmx) expzmid = exp(-.5*dz) expzmid_inv = 1./expzmid do i=1,3 rmassinv(i) = 1./rmass(i) enddo js=-jmaxh do j=1, jmax phi=(j+js-.5)*dphi cs(j)=cos(phi) sn(j)=sin(phi) tn(j)=tan(phi) cor(j)=2.*omega*sn(j) enddo ! ! cs at 0, -1, zjmx+1, and zjmx+2 replace the old cssp and csnp: cs(-1) = -cs(2) cs(0) = -cs(1) cs(zjmx+1) = -cs(zjmx) cs(zjmx+2) = -cs(zjmx-1) dt = float(step) ! was C(4) dtx2 = 2.*dt ! was C(6) dtx2inv = 1./dtx2 ! was C(7) freq_3m3 = 2.*pi/(49.7789*60.*60.) ! was C(21) freq_semidi = 4.*pi/(24.*60.*60.) ! was C(23) freq_ann = freq_semidi/(2.*365.25) ! was C(25) ! ! Set gravity according to lower boundary: if (zsb==-17.) then ! time-gcm grav = 945. elseif (zsb==-7.) then ! tiegcm grav = 870. else write(6,"(/,'>>> WARNING: do not know how to assign gravity', | ' constant with zsb=',f10.2)") zsb grav = 945. write(6,"(' Default to grav=',f10.2,/)") grav endif dzgrav = grav/gask ! C(65) ! ! Report to stdout: write(6,"(/,'Model name = ',a)") tgcm_name write(6,"( 'Model version = ',a)") tgcm_version write(6,"(/,'Set constants:')") write(6,"(' zjmx=',i3,' zimx=',i3,' zkmx=',i3)") zjmx,zimx,zkmx write(6,"(' zst=',f8.2,' zsb=',f8.2,' dz=',f5.2)") zst,zsb,dz write(6,"(' zfldx=',i5,' zphyx=',i5)") zfldx,zphyx write(6,"(' dt = ',f8.2,' secs')") dt write(6,"(' grav = ',f10.2)") grav write(6,"(' freq_3m3 = ',e10.4,' freq_semidi=',e10.4, | ' freq_ann=',e12.4)") freq_3m3,freq_semidi,freq_ann write(6,"(' dipmin = ',f8.3)") dipmin ! end subroutine init_cons !----------------------------------------------------------------------- subroutine set_consdyn ! ! Define constants for dynamo (consdyn.h): ! implicit none #include "params.h" #include "consdyn.h" real :: tabl2,tanth0,tanths,theta0,hamh0 COMMON/BLNK/TABL2(91,3:5),TANTH0(JMAXM),TANTHS(JMAXM), 1 THETA0(JMAXM),HAMH0(JMAXM) ! ! constants in consdyn.h common: ! DATA RE/6.378165E8/ ! DATA H0/9.0E6/ ! DATA H00/9.7E6/ ! DATA HS/1.3E7/ ! DATA UNIT/IMAXM*1.0/ ! DATA R1/1.06E7/ ! ! Local: integer :: i,j,n real :: e,alfa,dtheta,tanths2 ! ! constants /consdyn/ in consdyn.h: ! (should be parameter constants in this module, ! then use-associated to other routines, but ! note re above is different this one) ! RE = 6.378165E8 H0 = 9.0E6 H00 = 9.7E6 HS = 1.3E7 R1 = 1.06E7 UNIT(:)=1.0 ! E=1.E-6 ALFA = 1.668 R0 = RE+H0 R00 = RE+H00 RS = RE+HS ! PI = 4.*ATAN(1.) ! pi is in this module DTR = PI/180. RTD = 180./PI DLATG = PI/FLOAT(JMAXG) DLONG = 2.*PI/FLOAT(IMAXG) DLATM = PI/FLOAT(JMAXM-1) DLONM = 2.*PI/FLOAT(IMAXM) C **** C **** FILL ARRAY YLATG C **** DO 1 J = 1,JMAXG YLATG(J)=-.5*(PI-DLATG)+FLOAT(J-1)*DLATG 1 CONTINUE YLATG(0) = -PI/2.+E YLATG(JMAXGP) = PI/2.-E C **** C **** FILL YLONG C **** DO 2 I = 1,IMAXGP YLONG(I) = -PI+FLOAT(I-1)*DLONG 2 CONTINUE C **** C **** FILL ARRAY YLATM (EQUALLY SPACED IN THETAO BUT HOLDS C **** CORRESPONDING VALUE OF THETAS) C **** DO 3 J = 1,JMAXM C **** THETA0 = EQUALLY SPACED GRID VALUES THETA0(J) = -PI/2.+FLOAT(J-1)*DLATM 3 CONTINUE DO 4 J = 2,JMAXM-1 C **** TANTH0 = ABS(TAN(THETA0)) TANTH0(J) = ABS(TAN(THETA0(J))) C **** HAMH0 = HA-H0 HAMH0(J) = R1*TANTH0(J)+R0*TANTH0(J)**(2.+2.*ALFA)/ 1 (1.+TANTH0(J)**2)**ALFA C **** TANTHS = ABS(TAN(THETAS)) TANTHS(J) = SQRT(HAMH0(J)/R0) C **** YLATM = TANTHS YLATM(J) = SIGN(ATAN(TANTHS(J)),THETA0(J)) C **** RCOS0S = COS(THETA0)/COS(THETAS) RCOS0S(J) = SQRT((1.+TANTHS(J)**2)/(1.+TANTH0(J)**2)) C **** DT0DTS = D(THETA0)/D(THETAS) c calculate DT1DTS C DT1DTS = DT0DTS divided by abs[sin(Im)]. C Remains finite and nonzero at magnetic equator tanths2 = TANTHS(J)**2 DT1DTS(J) = (R0*sqrt(1.+4.*tanths2)*(1.+tanths2))/ 1 (R1*(1.+TANTH0(J)**2)+2.*R0*TANTH0(J)**(2.*ALFA+1.)* 2 (1.+ALFA+TANTH0(J)**2)/(1.+TANTH0(J)**2)**ALFA) DT0DTS(J) = DT1DTS(J)*2.*TANTHS(J)/sqrt(1.+4.*tanths2) 4 CONTINUE C **** C **** NOW DO POLES C **** YLATM(1) = THETA0(1) YLATM(JMAXM) = THETA0(JMAXM) RCOS0S(1) = 1. RCOS0S(JMAXM) = 1. DT0DTS(1) = 1. DT0DTS(JMAXM) = 1. DT1DTS(1) = 1. DT1DTS(JMAXM) = 1. C **** C **** FILL YLONM C **** DO 5 I = 1,IMAXMP YLONM(I) = -PI+FLOAT(I-1)*DLONM 5 CONTINUE DTHETA = PI/(2.*90.) TABLE(1,1) = 0. TABLE(1,2) = 0. DO 6 I = 2,91 TABLE(I,1) = TABLE(I-1,1)+DTHETA 6 CONTINUE DO 7 I = 2,90 TABL2(I,4) = TAN(TABLE(I,1)) TABLE(I,2) = TABLE(I,1) 7 CONTINUE DO 8 N = 1,7 DO 9 I = 2,90 TABL2(I,3) = TABLE(I,2) TABLE(I,2) = TAN(TABL2(I,3)) TABL2(I,5) = SQRT(R1/R0*TABLE(I,2)+TABLE(I,2)** 1 (2.*(1.+ALFA))/(1.+TABLE(I,2)**2)**ALFA) TABLE(I,2) = TABL2(I,3)-(TABL2(I,5)-TABL2(I,4))*2.*TABL2(I,5)/ 2 (R1/R0*(1.+TABLE(I,2)**2)+2.*TABLE(I,2)**(2.*ALFA+1.)* 3 (1.+ALFA+TABLE(I,2)**2)/(1.+TABLE(I,2)**2)**ALFA) 9 CONTINUE 8 CONTINUE RETURN end subroutine set_consdyn end module cons_module