module proc_lat implicit none contains !----------------------------------------------------------------------- subroutine proclat(flat,fproc,f,nf,h,imx,kmx,lat,modelhts,iprint) ! ! Processing for "new" history file format (11/05 and later for tiegcm): ! use hist,only: history ! use only history type-def use fields,only: field ! use only field type-def use netcdf_module,only: tlbc,ulbc,vlbc use proc,only: nlon,gcmlat,spval ! ! Args: type(history),intent(in) :: h ! the current history integer,intent(in) :: nf ! number of field structs type(field),intent(inout) :: f(nf) ! field structs 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) ! ! Local: integer :: i,ix,ixproc,ixf,k,ixo1,ixo2,ixt,ixz,ixco,ixn2 real,dimension(imx,kmx) :: raw_tn, raw_o1, raw_o2, raw_z, raw_co, | raw_n2 real,parameter :: eps=1.e-12 real,parameter :: boltz=1.38e-16 ! 2/14/08 btf: change lower bound of venus Z as per Bougher: ! real,parameter :: cmbot_vtgcm = 94.0*1.e+5 ! bottom Z for vtgcm (cm) real,parameter :: cmbot_vtgcm = 69.0e+05 ! bottom Z for vtgcm (cm) ! real :: pzps(kmx) ! scale heights ! ! Externals: integer,external :: ixfindc ! ! 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) if (h%isvtgcm) then ixco = ixfindc(h%fnames,h%nflds,'CO ') if (ixco==0) then write(6,"('>>> proclat: need CO in flat for vtgcm')") stop 'CO' endif raw_co(:,:) = flat(:,:,ixco) ! write(6,"('proclat vtgcm: lat=',i3,' ixco=',i3, ! | ' co min,max=',2e12.4)") lat,ixco, ! | minval(raw_co),maxval(raw_co) ! ixn2 = ixfindc(h%fnames,h%nflds,'N2 ') if (ixn2==0) then write(6,"('>>> proclat: need N2 in flat for vtgcm')") stop 'N2' endif raw_n2(:,:) = flat(:,:,ixn2) ! write(6,"('proclat vtgcm: lat=',i3,' ixn2=',i3, ! | ' n2 min,max=',2e12.4)") lat,ixn2, ! | minval(raw_n2),maxval(raw_n2) endif ! fields_loop: do ix=1,h%nflds ixproc = ixfindc(fproc,h%nflds,h%fnames(ix)) ixf = ixfindc(f%fname8,nf,h%fnames(ix)) ! write(6,"('proclat: field ',a,' ix=',i3,' ixproc=',i3,' ixf=', ! | i3)") h%fnames(ix),ix,ixproc,ixf if (ixproc == 0) cycle fields_loop ! ! Select on field to be processed: ! select case (h%fnames(ix)) ! ! For TN,UN,VN, the default behavior is to shift to interfaces, and use ! TLBC,ULBC,VLBC for the bottom interface boundary. If this is overridden ! by user request (zptype_req), then leave them at midpoints. ! (assume values at kmx are missing, do not change them) ! If switching zptype is requested for other fields (including derived), ! it will be done later by sub set_zptype (getflds.F). ! case ('TN ') if (ixf==0) then call mkinterfaces(flat(:,:,ix),tlbc(:,lat),imx,kmx) elseif (trim(f(ixf)%zptype_req)=='INTERFACES') then call mkinterfaces(flat(:,:,ix),tlbc(:,lat),imx,kmx) else ! leave at midpoints f(ixf)%zptype = 'MIDPOINTS' f(ixf)%lev = h%lev ! if (lat==1) write(6,"('proclat leaving TN at midpoints')") endif ! ! For UN, use ULBC for bottom boundary, and shift to interfaces: ! (assume values at kmx are missing, do not change them) case ('UN ') if (ixf==0) then call mkinterfaces(flat(:,:,ix),ulbc(:,lat),imx,kmx) flat(1:nlon,1:kmx-1,ix) = flat(1:nlon,1:kmx-1,ix)*.01 ! cm to m elseif (trim(f(ixf)%zptype_req)=='INTERFACES') then call mkinterfaces(flat(:,:,ix),ulbc(:,lat),imx,kmx) flat(1:nlon,1:kmx-1,ix) = flat(1:nlon,1:kmx-1,ix)*.01 ! cm to m else ! leave at midpoints f(ixf)%zptype = 'MIDPOINTS' f(ixf)%lev = h%lev flat(:,1:kmx-1,ix) = flat(:,1:kmx-1,ix)*.01 ! cm to m endif if (.not.all(flat(:,kmx,ix)==0.).and. | .not.all(flat(:,kmx,ix)==spval)) | flat(:,kmx,ix) = flat(:,kmx,ix)*.01 ! ! For VN, use VLBC for bottom boundary, and shift to interfaces: ! (assume values at kmx are missing, do not change them) case ('VN ') if (ixf==0) then call mkinterfaces(flat(:,:,ix),vlbc(:,lat),imx,kmx) flat(1:nlon,1:kmx-1,ix) = flat(1:nlon,1:kmx-1,ix)*.01 ! cm to m elseif (trim(f(ixf)%zptype_req)=='INTERFACES') then call mkinterfaces(flat(:,:,ix),vlbc(:,lat),imx,kmx) flat(1:nlon,1:kmx-1,ix) = flat(1:nlon,1:kmx-1,ix)*.01 ! cm to m else ! leave at midpoints f(ixf)%zptype = 'MIDPOINTS' f(ixf)%lev = h%lev flat(:,1:kmx-1,ix) = flat(:,1:kmx-1,ix)*.01 ! cm to m ! if (lat==1) write(6,"('proclat leaving VN at midpoints')") endif if (.not.all(flat(:,kmx,ix)==0.).and. | .not.all(flat(:,kmx,ix)==spval)) | flat(:,kmx,ix) = flat(:,kmx,ix)*.01 ! ! Recalculate Z using varying gravity if requested. Convert from cm to km. case ('Z ') if (modelhts > 0) then ! use heights on history flat(:,:,ix) = flat(:,:,ix) * 1.e-5 ! cm to km else if (h%isvtgcm) then call vcalchts(cmbot_vtgcm) else call calchts ! recalculate using varying gravity endif endif ! ! ZG is geopotential with varying gravity (only on new histories): case ('ZG ') flat(:,:,ix) = flat(:,:,ix) * 1.e-5 ! ! CO: case ('CO ') if (h%isvtgcm) then call mvtgcm_lbo1co('CO',h%version) endif ! ! O1: case ('O1 ') if (h%isvtgcm) then call mvtgcm_lbo1co('O1',h%version) endif case default end select ! write(6,"('proclat: field ',a,' zptype=',a,' lev=',/, ! | (9f8.3))") f(ixf)%fname8,f(ixf)%zptype,f(ixf)%lev enddo fields_loop return ! return from proclat ! ! Subroutines internal to proclat: contains !------------------------------------------------------------------- subroutine mkinterfaces(f,lbc,imx,kmx) use proc,only: nlon ! ! Args: integer,intent(in) :: imx,kmx real,intent(inout) :: f(imx,kmx) real,intent(in) :: lbc(imx) ! ! Local: integer :: k do k=kmx-1,2,-1 f(:,k) = 0.5*(f(:,k)+f(:,k-1)) enddo f(1:nlon-1,1) = lbc(:) f(nlon,1) = f(1,1) ! periodic point end subroutine mkinterfaces !------------------------------------------------------------------- subroutine calchts ! ! 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 !------------------------------------------------------------------- subroutine glatf(lat,gv,reff) ! ! CALCULATE LATITUDE VARIABLE GRAVITY (GV) AND EFFECTIVE ! RADIUS (REFF) ! 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 vcalchts(zcmbot) ! ! Calculate heights as in calchts above, but for vtgcm: ! (Latest version : 07/10/00 : S. W. BOUGHER) ! (Latest version : 10/12/00 : B. FOSTER & S. W. BOUGHER) ! 2/14/08 btf: changed cmbot_vtgcm (zcmbot) and g0 as per 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 ! ! 10/27/06 btf: parameters for venus vtgcm: ! 2/14/08 btf: change g0 as per Bougher: ! real,parameter :: g0=887.0, r=6.052e+8, dz=0.5 ! 1/7/11 btf: Calculate dz rather than hardwiring, so this will ! work for any vertical resolution. real,parameter :: g0=870.0, r=6.052e+8 real :: g(kmx),fco2(kmx),xmas(kmx),dz 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 vcalchts !------------------------------------------------------------------- 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 ! 2/14/08 btf: change co1,cco according to correspondence from Bougher ! co1 = 2.0e-6 ! cco = 1.5e-3 co1 = 9.35e-10 cco = 4.56e-05 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) ! write(6,"(/,'mvtgcm_lbo1co ',a,' flat before pso=',/, ! | (6e12.4))") name,flat(:,1,ix) ! write(6,"('mvtgcm_lbo1co ',a,': fb=',e12.4,' pso=',/, ! | (6e12.4))") name,fb,pso flat(:,1,ix) = 0.5*(flat(:,1,ix)+pso(:)) ! write(6,"('mvtgcm_lbo1co ',a,': final flat lbc=',/,(6e12.4))") ! | name,flat(:,1,ix) end subroutine mvtgcm_lbo1co !------------------------------------------------------------------- end subroutine proclat end module proc_lat