module proc_lat_old implicit none contains !----------------------------------------------------------------------- subroutine proclat_old(flat,fproc,h,imx,kmx,lat,modelhts, | iprint) ! ! Redefine flat(h%nlon-1,h%nzp,h%nflds) from midpoint to interface levels ! for requested fields fproc. ! Change units of non-density species where necessary. ! Calculate heights if modelhts=0 ! ! flat has been read from history for the current latitude index lat. ! (imx==nlon===h%nlon-1, kmx==npress==h%nzp) ! Field names in lat slice are in history_type h (h%fnames). ! lat is the current latitude index. ! This routine does the following as needed: ! -- interpolation/extrapolation from midpoint to interface levels ! -- unit conversions (except densities) ! use hist,only: history ! use only history type-def use fields,only: flds_known,nfknown use input,only: zkmbot_mgcm, tbot_mgcm use proc,only: gcmlat ! ! Args: type(history),intent(in) :: h ! the current history integer,intent(in) :: imx,kmx ! number of lons and zp levels real,intent(inout) :: flat(imx,kmx,h%nflds) ! the current lat slice character(len=8),intent(in) :: fproc(h%nflds) ! fields to be processed integer,intent(in) :: lat ! current lat index integer,intent(in) :: modelhts ! heights flag (from input) integer,intent(in) :: iprint ! print flag (currently levels 0 or 1) ! ! Locals: real,parameter :: boltz=1.38e-16, eps=1.e-12 integer :: k,i,ix, ! loop indices + ixt,ixo2,ixo1,ixz,ixco, ! indices to selected fields in flat | ixh,ixhe,ixn2, + ixknown, ! index to flds_known if current field ! is known to the processor + ixproc, ! > 0 if field is to be processed + ier real :: pzps(kmx) ! scale heights real :: raw_tn(imx,kmx),raw_o1(imx,kmx),raw_o2(imx,kmx), + raw_z(imx,kmx), ! unprocessed t,o1,o2,z | raw_h(imx,kmx), raw_he(imx,kmx) real,allocatable :: raw_co(:,:) ! unprocessed co (needed for mtgcm) real,allocatable :: raw_n2(:,:) ! unprocessed n2 (needed for mtgcm) real :: fmin,fmax,cmbot_mtgcm ! ! Externals: integer,external :: ixfindc ! ! Check dimensions: ! if (imx /= h%nlon-1) then write(6,"('WARNING proclat: imx /= h%nlon-1: imx=',i3, + ' h%nlon=',i3)") imx,h%nlon stop 'proclat' endif if (kmx /= h%nzp) then write(6,"('WARNING proclat: kmx /= h%nzp: kmx=',i3, + ' h%nzp=',i3)") kmx,h%nzp stop 'proclat' endif ! ! Define needed field indices in flat: ! (these fields were confirmed on the history by verify_hist) ! ixo1 = ixfindc(h%fnames,h%nflds,'O1 ') where(flat(:,:,ixo1) < 0.) flat(:,:,ixo1) = eps raw_o1(:,:) = flat(:,:,ixo1) ixo2 = ixfindc(h%fnames,h%nflds,'O2 ') raw_o2(:,:) = flat(:,:,ixo2) ixt = ixfindc(h%fnames,h%nflds,'TN ') raw_tn(:,:) = flat(:,:,ixt) ixz = ixfindc(h%fnames,h%nflds,'Z ') raw_z(:,:) = flat(:,:,ixz) ! ! Save CO and N2 for mtgcm: if (h%ismtgcm) then ixco = ixfindc(h%fnames,h%nflds,'CO ') if (ixco==0) then write(6,"('>>> proclat: need CO in flat for mtgcm')") stop 'CO' endif allocate(raw_co(imx,kmx),stat=ier) if (ier /= 0) call allocerr(ier,"proclat allocate raw_co") raw_co(:,:) = flat(:,:,ixco) ! ixn2 = ixfindc(h%fnames,h%nflds,'N2 ') if (ixn2==0) then write(6,"('>>> proclat: need N2 in flat for mtgcm')") stop 'N2' endif allocate(raw_n2(imx,kmx),stat=ier) if (ier /= 0) call allocerr(ier,"proclat allocate raw_n2") raw_n2(:,:) = flat(:,:,ixn2) endif ! ismtgcm ! ! For jtgcm, sub jcalchts: ! 6/3/05: added conditionals on ixh,ixhe to satisfy latest Linux ! kernel at hao. ixh = ixfindc(h%fnames,h%nflds,'H ') if (ixh > 0) raw_h(:,:) = flat(:,:,ixh) ixhe = ixfindc(h%fnames,h%nflds,'HE ') if (ixhe > 0) raw_he(:,:) = flat(:,:,ixhe) ! ! Loop through fields, redefining flat(:,:,ix) if needed. ! (field is needed if h%fnames(ix) is in fproc list) ! fields_loop: do ix=1,h%nflds ixproc = ixfindc(fproc,h%nflds,h%fnames(ix)) if (ixproc == 0) then ! write(6,"('proclat: ix=',i3,' field h%fnames(ix)=',a, ! | ' apparently not needed (not in fproc)')") ix,h%fnames(ix) cycle fields_loop endif ixknown = ixfindc(flds_known%fname8,nfknown,h%fnames(ix)) ! ! Select on field to be processed: ! select case (h%fnames(ix)) ! ! TN,UN,VN: interpolate body to interface levels. Bottom level is stored ! in top slot, and kmx gets kmx-1. ! Note interp_body is called *before* setting bottom boundary from top slot. ! Interp_body sets interface levels from k==2 to k==kmx-1. Calling interp_body ! before setting lbc insures that interp_body uses the bottom midpoint-level ! at k==1 rather than the bottom boundary itself when setting interface level ! k==2. ! case ('TN ') if (.not.h%ismtgcm) then ! (including VTGCM) ! if (h%version(1:4) == 'tgcm') ! | flat(:,1,ix) = flat(:,kmx,ix) call interp_body(kmx-1) flat(:,1,ix) = flat(:,kmx,ix) else ! ! Lower boundary of tn for mtgcm is from NASA AMES GCM (see getmgcm): ! If tbot_mgcm(lat)==0., then assume coupled mtgcm/mgcm history, and ! use values from the history, otherwise tbot_mgcm was read from disk ! file (see getmgcm). ! call interp_body(kmx-1) if (tbot_mgcm(lat) /= 0.) then ! mtgcm standalone flat(:,1,ix) = tbot_mgcm(lat)+(flat(:,kmx,ix)-140.) endif endif flat(:,kmx,ix) = flat(:,kmx-1,ix) ! ! UN and VN: interpolate body and convert from cm/s to m/s: case ('UN ') call interp_body(kmx-1) flat(:,1,ix) = flat(:,kmx,ix) flat(:,kmx,ix) = flat(:,kmx-1,ix) flat(:,:,ix) = flat(:,:,ix)*.01 case ('VN ') call interp_body(kmx-1) flat(:,1,ix) = flat(:,kmx,ix) flat(:,kmx,ix) = flat(:,kmx-1,ix) flat(:,:,ix) = flat(:,:,ix)*.01 ! ! O2: interpolate body (from kmx) and extrapolate bottom boundary: case ('O2 ') call interp_body(kmx) call extrap_bot ! ! OX (O+O3): case ('OX ') call interp_body(kmx) call extrap_bot ! ! N4S: case ('N4S ') call interp_body(kmx) call extrap_bot ! ! NOZ (NO+NO2): case ('NOZ ') call interp_body(kmx) call extrap_bot ! ! CO: case ('CO ') call interp_body(kmx) if (.not.h%ismtgcm.and.trim(h%version)/='VTGCM') then call extrap_bot else call mvtgcm_lbo1co('CO',h%version) endif ! ! CO2 (if mtgcm, CO2 is a derived field, see mkderived): case ('CO2 ') call interp_body(kmx) call extrap_bot ! ! H2O: case ('H2O ') call interp_body(kmx) call extrap_bot ! ! H2: case ('H2 ') call interp_body(kmx) call extrap_bot ! ! HOX (OH+HO2+H): case ('HOX ') call interp_body(kmx) call extrap_bot ! ! O+: log interpolate body (from kmx-1) and log extrapolate boundaries: case ('O+ ') flat(:,kmx,ix) = sqrt(flat(:,kmx-1,ix)**3 / + flat(:,kmx-2,ix)) do k=kmx-1,2,-1 flat(:,k,ix) = sqrt(flat(:,k,ix)*flat(:,k-1,ix)) enddo flat(:,1,ix) = flat(:,1,ix)**2 / flat(:,2,ix) ! ! OP: log interpolate body (from kmx-1) and log extrapolate boundaries: ! (OP is an alias for O+ on netcdf histories) case ('OP ') flat(:,kmx,ix) = sqrt(flat(:,kmx-1,ix)**3 / + flat(:,kmx-2,ix)) do k=kmx-1,2,-1 flat(:,k,ix) = sqrt(flat(:,k,ix)*flat(:,k-1,ix)) enddo flat(:,1,ix) = flat(:,1,ix)**2 / flat(:,2,ix) ! ! CH4: case ('CH4 ') call interp_body(kmx) call extrap_bot ! ! AR: case ('AR ') call interp_body(kmx) call extrap_bot ! ! HE: case ('HE ') call interp_body(kmx) if (.not.h%isjtgcm) then call extrap_bot else ! is jtgcm call jtgcm_lbheh('HE') endif ! ! NAT: case ('NAT ') call interp_body(kmx) call extrap_bot ! ! O21D: case ('O21D ') call extrap_top call interp_body(kmx-1) call extrap_bot ! ! NO2: case ('NO2 ') call extrap_top call interp_body(kmx-1) call extrap_bot ! ! NO: case ('NO ') call extrap_top call interp_body(kmx-1) call extrap_bot ! ! O3: case ('O3 ') call extrap_top call interp_body(kmx-1) call extrap_bot ! ! O1: ! Note: may not need to extrapolate top boundary -- model ! does calculate top (see tgcm21 comp.f.259-264) case ('O1 ') if (.not.h%ismtgcm.and.trim(h%version)/='VTGCM') then call extrap_top ! maybe not necessary? call interp_body(kmx-1) call extrap_bot else call interp_body(kmx) call mvtgcm_lbo1co('O1',h%version) endif ! ! N2 (mtgcm only): case ('N2 ') if (.not.h%ismtgcm) then call interp_body(kmx) call extrap_bot endif ! ! OH: case ('OH ') call extrap_top call interp_body(kmx-1) call extrap_bot ! ! HO2: case ('HO2 ') call extrap_top call interp_body(kmx-1) call extrap_bot ! ! H: case ('H ') if (.not.h%isjtgcm) then call extrap_top call interp_body(kmx-1) call extrap_bot else call interp_body(kmx) call jtgcm_lbheh('H') endif ! ! N2D: case ('N2D ') call extrap_top call interp_body(kmx-1) call extrap_bot ! ! TI: case ('TI ') call extrap_top call interp_body(kmx-1) call extrap_bot ! ! TE: case ('TE ') call extrap_top call interp_body(kmx-1) call extrap_bot ! ! NE (no change): case ('NE ') ! ! O2+: log interpolate body (from kmx-1) and log extrapolate boundaries: ! (check for o2+ <= eps at boundaries before extrapolating) case ('O2+ ') where(flat(:,:,ix) < 0.) flat(:,:,ix) = eps flat(:,kmx,ix) = sqrt(flat(:,kmx-1,ix)**3 / flat(:,kmx-2,ix)) do k=kmx-1,2,-1 flat(:,k,ix) = sqrt(flat(:,k,ix)*flat(:,k-1,ix)) enddo flat(:,1,ix) = flat(:,1,ix)**2 / flat(:,2,ix) ! ! (O2P is an alias for O+ on netcdf histories) case ('O2P ') where(flat(:,:,ix) < 0.) flat(:,:,ix) = eps flat(:,kmx,ix) = sqrt(flat(:,kmx-1,ix)**3 / flat(:,kmx-2,ix)) do k=kmx-1,2,-1 flat(:,k,ix) = sqrt(flat(:,k,ix)*flat(:,k-1,ix)) enddo flat(:,1,ix) = flat(:,1,ix)**2 / flat(:,2,ix) ! ! W = vertical velocity: case ('W ') do i=1,imx do k=2,kmx-1 pzps(k) = raw_z(i,k+1)-raw_z(i,k-1) enddo pzps(1) = 4.*raw_z(i,2)-3.*raw_z(i,1)-raw_z(i,3) pzps(kmx) = 3.*raw_z(i,kmx)-4.*raw_z(i,kmx-1)+raw_z(i,kmx-2) flat(i,:,ix) = flat(i,:,ix)*pzps(:) enddo flat(imx,:,ix) = flat(1,:,ix) ! force periodic point flat(:,:,ix) = flat(:,:,ix) * .01 ! cm/s to m/s ! ! Z = heights ! If modelhts==0, calculate heights from mean mass and ht-varying gravity ! If modelhts==1, get model heights from history (may need to add zkmbot ! for old style histories -- check lower boundary) ! (zkmbot_mgcm(nlat) is bottom boundary of Z in mtgcm, read from nasa ! ames data file. See sub getmgcm in gethist.F) ! case ('Z ') ! write(6,"('proclat Z: modelhts=',i3,' h%version=',a, ! | ' zkmbot_mgcm=',/,(6e12.4))") modelhts,h%version,zkmbot_mgcm if (modelhts > 0) then ! get heights from history ! ! If standalone mtgcm run, zkmbot_mgcm(lat) was read from disk file, ! otherwise (coupled mtgcm/mgcm), it is all zeroes: ! if (h%ismtgcm) then flat(:,:,ix) = flat(:,:,ix)*1.e-5 + zkmbot_mgcm(lat) elseif (trim(h%version)=='VTGCM') then flat(:,:,ix) = flat(:,:,ix)*1.e-5 + 94. elseif (h%isjtgcm) then flat(:,:,ix) = flat(:,:,ix)*1.e-5 + 250. else fmax = maxval(flat(:,1,ix)) ! max at bottom if (fmax*1.e-5 > -1. .and. fmax*1.e-5 < 1.) then flat(:,:,ix) = flat(:,:,ix) * 1.e-5 + 97. else flat(:,:,ix) = flat(:,:,ix) * 1.e-5 ! do k=1,kmx ! write(6,"('proclat_old Z: lat=',i3,' k=',i3, ! | ' flat(:,k,ix)=',/,(6e12.4))") lat,k, ! | flat(:,k,ix) ! enddo endif endif else ! calculate heights if (h%ismtgcm) then ! ! Average bottom boundary of Z for mtgcm will be used by mcalchts ! only if zkmbot_mgcm is all 0. This is the case for a coupled ! mtgcm/mgcm run (bottom boundaries were received directly from mgcm). ! ! For mtgcm/mgcm coupled history: ! zkmbot_mgcm(:)==0. (Z bottom boundary NOT read from disk file) ! Mcalchts uses zonal mean Z from history for bottom boundary. ! For mtgcm standalone history: ! zkmbot_mgcm(lat) bottom boundary was read from disk file (see getmgcm) ! Mcalchts uses zkmbot_mgcm(lat) ! if (zkmbot_mgcm(lat)==0.) then ! coupled mtgcm/mgcm cmbot_mtgcm = 0. do i=1,imx cmbot_mtgcm=cmbot_mtgcm+flat(i,1,ix) ! already in cm enddo cmbot_mtgcm=cmbot_mtgcm/float(imx) else ! mtgcm standalone cmbot_mtgcm = zkmbot_mgcm(lat)*1.e5 ! km->cm endif call mcalchts(cmbot_mtgcm) elseif (h%isjtgcm) then call jcalchts else call calchts_old endif endif ! ! POTEN = electric potential (no change): case ('POTEN ') ! ! ENO (vtgcm only): case ('ENO ') call interp_body(kmx) call extrap_bot ! ! jtgcm fields H3P, HP, H2V1-V4, H2P, MORPH. ! ! H3+ case ('H3P ') flat(:,1,ix) = flat(:,2,ix) call interp_body(kmx-1) flat(:,kmx,ix) = flat(:,kmx-1,ix) ! ! H+ case ('HP ') flat(:,1,ix) = flat(:,2,ix) call interp_body(kmx-1) flat(:,kmx,ix) = flat(:,kmx-1,ix) ! ! H2V1-V4 case ('H2V1 ') flat(:,1,ix) = flat(:,2,ix) call interp_body(kmx-1) flat(:,kmx,ix) = flat(:,kmx-1,ix) case ('H2V2 ') flat(:,1,ix) = flat(:,2,ix) call interp_body(kmx-1) flat(:,kmx,ix) = flat(:,kmx-1,ix) case ('H2V3 ') flat(:,1,ix) = flat(:,2,ix) call interp_body(kmx-1) flat(:,kmx,ix) = flat(:,kmx-1,ix) case ('H2V4 ') flat(:,1,ix) = flat(:,2,ix) call interp_body(kmx-1) flat(:,kmx,ix) = flat(:,kmx-1,ix) ! ! H2+ case ('H2P ') flat(:,1,ix) = flat(:,2,ix) call interp_body(kmx-1) flat(:,kmx,ix) = flat(:,kmx-1,ix) ! ! MORPH case ('MORPH ') flat(:,1,ix) = flat(:,2,ix) call interp_body(kmx-1) flat(:,kmx,ix) = flat(:,kmx-1,ix) ! ! Unknown field: case default if (iprint > 0) + write(6,"('>>> proclat: do not recognize field ',a)") + h%fnames(ix) end select enddo fields_loop if (h%ismtgcm) deallocate(raw_co) return ! return from proclat ! ! Subroutines internal to proclat: contains !------------------------------------------------------------------- subroutine interp_body(itop) ! ! Interpolate body from k=itop down to k=2 (densities should not go < 0): ! integer,intent(in) :: itop integer :: k do k=itop,2,-1 flat(:,k,ix) = 0.5*(flat(:,k,ix)+flat(:,k-1,ix)) if (ixknown > 0 .and. + flds_known(ixknown)%type == 'DENSITY ' .and. + minval(flat(:,k,ix)) < 0. .and. iprint > 0) + write(6,"('WARNING interp_body: ',a,' < 0.: lat=',i2, + ' min,max of lon values at k=',i2,': ',2e12.4)") + flds_known(ixknown)%fname8,lat,k,minval(flat(:,k,ix)), + maxval(flat(:,k,ix)) enddo return end subroutine interp_body !------------------------------------------------------------------- subroutine extrap_top ! ! Extrapolate to get top boundary (insure positive densities): ! (o3 often goes < 0 here) ! real :: av integer :: i flat(:,kmx,ix) = 1.5*flat(:,kmx-1,ix)-0.5*flat(:,kmx-2,ix) if (ixknown > 0 .and. + flds_known(ixknown)%type == 'DENSITY ') then if (minval(flat(:,kmx,ix)) < 0.) then ! av = sum(abs(flat(:,kmx,ix))) / float(imx-1) ! where(flat(:,kmx,ix) < 0.) flat(:,kmx,ix) = av do i=1,imx if (flat(i,kmx,ix) < 0.) flat(i,kmx,ix) = flat(i,kmx-1,ix) enddo if (iprint > 0) ! + write(6,"('extrap_top: lat=',i2,' field ',a, ! + ' min < 0. (ave=',e12.4,')')") ! + lat,flds_known(ixknown)%fname8,av + write(6,"('extrap_top: lat=',i2,' field ',a, + ' min < 0.')") lat,flds_known(ixknown)%fname8 endif endif end subroutine extrap_top !------------------------------------------------------------------- subroutine extrap_bot ! ! Extrapolate to get bottom boundary (insure positive densities): ! (o1 often < 0. here) ! ! In the following comments, primed numbers (e.g., 2') mean the ! midpoint level, as received from the history (original flat). ! Extrapolation of bottom boundary requires use of midpoint-levels ! 1 and 2 (f(1)=1.5*f(1')-.5f(2')), but midpoint-level 2 (2') ! was overwritten by interp_body (called before this routine). ! Thus, we retrieve midpoint-level 2 as follows: ! f(2) = 0.5*(f(2')+f(1') ==> f(2') = 2.*f(2)-f(1') ! This is then substituted for f(2') in the extrapolation to ! give f(1) = 2.*f(1')-f(2) ! real :: av integer :: i ! flat(:,1,ix) = 1.5*flat(:,1,ix)-0.5*flat(:,2,ix) flat(:,1,ix) = 2.*flat(:,1,ix)-flat(:,2,ix) if (ixknown > 0 .and. + flds_known(ixknown)%type == 'DENSITY ') then if (minval(flat(:,1,ix)) < 0.) then ! av = sum(abs(flat(:,1,ix))) / float(imx-1) ! where(flat(:,1,ix) < 0.) flat(:,1,ix)=av do i=1,imx if (flat(i,1,ix) < 0.) flat(i,1,ix) = flat(i,2,ix) enddo if (iprint > 0) ! + write(6,"('extrap_bot: lat=',i2,' field ',a, ! + ' min < 0. (ave=',e12.4,')')") ! + lat,flds_known(ixknown)%fname8,av + write(6,"('extrap_bot: lat=',i2,' field ',a, + ' min < 0.')") lat,flds_known(ixknown)%fname8 endif endif end subroutine extrap_bot !------------------------------------------------------------------- subroutine calchts_old ! ! Calculate heights from mean mass, with gravity varying w/ height: ! 5/25/05 btf: new sub calchts from liying (lqian), to include ! latitudinally varying gravity (sub glatf). This ! is included in tgcmproc1.7. ! real :: g0,r0 real :: g(kmx),xmas(kmx),dz real :: ftn(kmx),fo2(kmx),fo1(kmx),fn2(kmx) integer :: i dz=(h%zpt-h%zpb)/float(kmx-1) ! from cm to km ! do i=1,imx ! ! Process t,o2,o1,n2 and store in locals ftn,fo2,fo1,fn2: ! ftn,fo1,fo2,fn2 are values at model midpoint ftn(:)=raw_tn(i,:) fo1(:)=raw_o1(i,:) fo2(:)=raw_o2(i,:) fn2(:) = max(.00001,(1.-fo2(:)-fo1(:))) ! ! Get mean mass and calculate new heights: ! Note bottom boundary has already been read from the history. ! xmas and g are values at model midpoint as well xmas(:) = 1./(fo1(:)/16.+fo2(:)/32.+fn2(:)/28.)*1.66e-24 C calculate latitude dependent g0 (at Earth surface) and Earth radius r0 call glatf(gcmlat(lat),g0,r0) g(1)=g0*(r0/(r0+0.5*(flat(i,1,ix)+flat(i,2,ix))))**2 do k=2,kmx flat(i,k,ix) = flat(i,k-1,ix) + boltz*dz* + ftn(k-1) / (xmas(k-1)*g(k-1)) ! ! 3/22/06 btf: add conditional on k to avoid flat overflow at k+1 if (k < kmx) | g(k)=g0*(r0/(r0+0.5*(flat(i,k,ix)+flat(i,k+1,ix))))**2 enddo flat(i,:,ix) = flat(i,:,ix) * 1.e-5 ! cm to km enddo return end subroutine calchts_old !------------------------------------------------------------------- SUBROUTINE GLATF(LAT,GV,REFF) C CALCULATE LATITUDE VARIABLE GRAVITY (GV) AND EFFECTIVE C RADIUS (REFF) C Earth radius in cm real,intent(in) :: lat real,intent(out) :: gv,reff real :: dgtr=1.74533E-2 real :: c2 C2 = COS(2.*DGTR*LAT) GV = 980.616*(1.-.0026373*C2) REFF = 2.*GV/(3.085462E-6 + 2.27E-9*C2) RETURN end subroutine glatf !------------------------------------------------------------------- ! subroutine mcalchts(zcmbot) ! ! Calculate heights as in calchts above, but for mtgcm: ! ! real,intent(in) :: zcmbot ! bottom boundary in cm ! real,parameter :: g0=371.1, r=3388.25e+5 ! real :: dz,g(kmx),fn2(kmx),fco2(kmx),xmas(kmx) ! integer :: i,k ! ! dz=(h%zpt-h%zpb)/float(kmx-1) ! do i=1,imx ! write(6,"('raw_o2=',/,(6e12.4))") raw_o2(i,:) ! write(6,"('raw_o1=',/,(6e12.4))") raw_o1(i,:) ! write(6,"('raw_co=',/,(6e12.4))") raw_co(i,:) ! fn2(:) = max(.00001,(1.-raw_o2(i,:)-raw_o1(i,:))) ! fco2(:) = max(.00001,(1.-raw_o1(i,:)-raw_co(i,:)-fn2(:))) ! g(1) = g0 ! xmas(:) = 1./(raw_o1(i,:)/16.+(raw_co(i,:)+fn2(:))/28.+ ! + fco2(:)/44.) * 1.66e-24 ! flat(i,1,ix) = zcmbot ! do k=2,kmx ! flat(i,k,ix) = flat(i,k-1,ix) + boltz*dz* ! + raw_tn(i,k-1) / (xmas(k-1)*g(k-1)) ! g(k) = g0*(r/(r+flat(i,k,ix)))**2 ! enddo ! flat(i,:,ix) = flat(i,:,ix) * 1.e-5 ! cm to km ! enddo ! end subroutine mcalchts !------------------------------------------------------------------- subroutine mcalchts(zcmbot) ! ! Calculate heights as in calchts above, but for mtgcm: ! (Latest version : 07/10/00 : S. W. BOUGHER) ! (Latest version : 10/12/00 : B. FOSTER & S. W. BOUGHER) ! (Always Using cm units for heights) ! real,intent(in) :: zcmbot ! bottom boundary in cm ! real,parameter :: g0=372.65, r=3388.25e+5 ! ------ MARS BOOK PARAMETERS (1992) MEAN VALUES ------------------- ! real,parameter :: g0=371.1, r=3388.25e+5, dz=0.5 ! ------ MOLA PARAMETERS (1999) MEAN VALUES ------------------------ ! real,parameter :: g0=371.1, r=3389.51e+5, dz=0.5 real,parameter :: g0=372.5, r=3389.51e+5, dz=0.5 ! real :: g(kmx),fco2(kmx),xmas(kmx) integer :: i,k ! dz=(h%zpt-h%zpb)/float(kmx-1) do i=1,imx fco2(:) = max(.00001,(1.-raw_o1(i,:)-raw_co(i,:)- + raw_n2(i,:))) xmas(:) = 1./(raw_o1(i,:)/16.+(raw_co(i,:)+raw_n2(i,:))/28.+ + fco2(:)/44.) * 1.66e-24 ! ----------------------------------------------- ! flat(i,1,ix) = zkmbot_mgcm(lat)*1.e5 ! Old Method flat(i,1,ix) = zcmbot ! ----------------------------------------------- g(1) = g0*(r/(r+flat(i,1,ix)))**2.0 do k=2,kmx flat(i,k,ix) = flat(i,k-1,ix) + boltz*dz* + raw_tn(i,k-1) / (xmas(k-1)*g(k-1)) g(k) = g0*(r/(r+flat(i,k,ix)))**2 enddo flat(i,:,ix) = flat(i,:,ix) * 1.e-5 ! cm to km enddo end subroutine mcalchts !------------------------------------------------------------------- subroutine mvtgcm_lbo1co(name,version) ! ! Set lower boundary of O1 and CO for mtgcm or vtgcm: ! character(len=*),intent(in) :: name,version real :: fb,pso(imx),co1,cco ! if (h%ismtgcm) then co1 = 5.0e-4 cco = 2.0e-3 elseif (trim(version)=='VTGCM') then co1 = 2.0e-6 cco = 1.5e-3 else write(6,"('>>> mvtgcm_lbo1co: unrecognized version=',a)")version write(6,"(' version must be MTGCM or VTGCM')") return endif if (trim(name)/='O1'.and.trim(name)/='CO') then write(6,"('>>> mvtgcm_lbo1co: unknown name=',a)") name return endif if (trim(name)=='O1') then fb = 2.*co1*(16./44.) else ! CO fb = 2.*cco*(28./44.) endif pso(:) = fb - flat(:,1,ix) flat(:,1,ix) = 0.5*(flat(:,1,ix)+pso(:)) end subroutine mvtgcm_lbo1co !------------------------------------------------------------------- subroutine jcalchts ! ! Calculate heights for jupiter gcm (not using geopotential from ! the history). Internal to sub proclat. ! ! Local: real :: | g(kmx), ! gravity accel | ftn(kmx), ! neutral temperature | barm(kmx), ! mean molecular weight | fhe(kmx), ! HE MMR from history | fh(kmx), ! H MMR from history | fh2(kmx), ! H2 MMR from history | hts(kmx) ! output heights real,parameter :: | go=2622., | ro=7.14e+9, | boltz=1.38e-16, ! boltman's constant | dz=0.5, ! delta Zp | htlb=250., ! height at bottom boundary | rgas=8.314E+07 ! gas constant ! do i=1,imx ftn(:) = raw_tn(i,:) fhe(:) = raw_he(i,:) fh(:) = raw_h(i,:) fh2(:) = 1.-fh(:)-fhe(:) ! ! Mean molecular weight from mixing ratios of He, H, H2: ! do k=1,kmx barm(k) = 1./(fhe(k)/4. + fh(k)/1. + fh2(k)/2.2) enddo ! ! Find heights (cm to km) ! hts(1) = htlb*1.e+5 g(1) = go*(ro/(ro+hts(1)))**2 do k=2,kmx hts(k) = hts(k-1) + rgas*dz*ftn(k-1) / (barm(k-1)*g(k-1)) g(k) = go*(ro/(ro+hts(k)))**2 enddo do k=1,kmx hts(k) = hts(k) * 1.e-5 enddo flat(i,:,ix) = hts(:) ! write(6,"('jcalchts: i=',i3,' hts=',/,(6e12.4))") i,hts enddo ! i=1,imx end subroutine jcalchts !------------------------------------------------------------------- subroutine jtgcm_lbheh(name) ! ! Set lower boundary of HE and H for jtgcm ! character(len=*),intent(in) :: name real :: fb,pso(imx),xhe,xh ! xhe = 1.35E-01 xh = 4.23E-08 ! if (trim(name)/='HE'.and.trim(name)/='H') then write(6,"('>>> jtgcm_lbheh: unknown name=',a)") name return endif if (trim(name)=='HE') then fb = 2.*xhe*(4./2.27) else ! H fb = 2.*xh*(1./2.27) endif pso(:) = fb - flat(:,1,ix) flat(:,1,ix) = 0.5*(flat(:,1,ix)+pso(:)) end subroutine jtgcm_lbheh end subroutine proclat_old !------------------------------------------------------------------- end module proc_lat_old