C +-======-+ C Copyright (c) 2003-2007 United States Government as represented by C the Admistrator of the National Aeronautics and Space Administration. C All Rights Reserved. C C THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, C REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN C COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS C REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). C THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN C INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR C REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, C DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED C HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE C RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. C C Government Agency: National Aeronautics and Space Administration C Government Agency Original Software Designation: GSC-15354-1 C Government Agency Original Software Title: GEOS-5 GCM Modeling Software C User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov C Government Agency Point of Contact for Original Software: C Dale Hithon, SRA Assistant, (301) 286-2691 C C +-======-+ program main implicit none c ********************************************************************** c ********************************************************************** c **** **** c **** Program to remap GEOS-5 FV & MOIST restarts in the vertical **** c **** **** c ********************************************************************** c ********************************************************************** character*256 dynrst, moistrst, bkgeta, topo integer headr1(6) integer headr2(5) integer nymd,nhms integer im,jm,lm_in,lm_out,nt,iostat,rc real undef, kappa, grav, getcon c restart variables and topography c -------------------------------- real*8, allocatable :: dp_in(:,:,:) real*8, allocatable :: u_in(:,:,:) real*8, allocatable :: v_in(:,:,:) real*8, allocatable :: thv_in(:,:,:) real*8, allocatable :: pk_in(:,:,:) real*8, allocatable :: ple_in(:,:,:) real*8, allocatable :: q_in(:,:,:,:) real*8, allocatable :: ps_in(:,:) real*8, allocatable :: ak_in(:) real*8, allocatable :: bk_in(:) real*8, allocatable :: phis(:,:) real*8, allocatable :: dp_out(:,:,:) real*8, allocatable :: u_out(:,:,:) real*8, allocatable :: v_out(:,:,:) real*8, allocatable :: thv_out(:,:,:) real*8, allocatable :: pk_out(:,:,:) real*8, allocatable :: pke_out(:,:,:) real*8, allocatable :: ple_out(:,:,:) real*8, allocatable :: q_out(:,:,:,:) real*8, allocatable :: ps_out(:,:) real*8, allocatable :: ak_out(:) real*8, allocatable :: bk_out(:) real*4, allocatable :: dum(:,:) real*4, allocatable :: zox(:,:) ! Zonal Mean odd-ox character*256, allocatable :: arg(:) character*8 date character*2 hour,clm character*3 cim,cjm integer n,nargs,iargc,i,j,L,nymd0,nhms0 C ********************************************************************** C **** Initialize Filenames **** C ********************************************************************** lm_out = -999 undef = 1.0e15 dynrst = 'fvcore_internal_restart' moistrst = 'moist_internal_restart' nargs = iargc() if(nargs == 0 ) call usage() allocate ( arg(nargs) ) do n=1,nargs call getarg(n,arg(n)) enddo do n=1,nargs if( trim(arg(n)).eq.'-h' ) call usage() if( trim(arg(n)).eq.'-help' ) call usage() if( trim(arg(n)).eq.'-H' ) call usage() if( trim(arg(n)).eq.'-Help' ) call usage() if( trim(arg(n)).eq.'-dynrst' ) then dynrst = trim(arg(n+1)) endif if( trim(arg(n)).eq.'-moistrst' ) then moistrst = trim(arg(n+1)) endif if( trim(arg(n)).eq.'-topo' ) then topo = trim(arg(n+1)) endif if( trim(arg(n)).eq.'-lm' ) read(arg(n+1),*) lm_out enddo if( lm_out.eq.-999 ) then print * print *, 'You must supply Output Vertical Resolution!' print * stop endif print * print *, ' dyn restart filename: ',trim(dynrst) print *, 'moist restart filename: ',trim(moistrst) print *, ' topo filename: ',trim(topo) C ********************************************************************** C **** Read dycore internal Restart **** C ********************************************************************** open (10,file=trim(dynrst),form='unformatted',access='sequential') read (10) headr1 read (10) headr2 nymd = headr1(1)*10000 + headr1(2)*100 + headr1(3) nhms = headr1(4)*10000 + headr1(5)*100 + headr1(6) im = headr2(1) jm = headr2(2) lm_in = headr2(3) print *, ' input resolution: ',im,jm,lm_in print *, ' date: ',nymd,nhms print * allocate ( u_in(im,jm,lm_in) ) allocate ( v_in(im,jm,lm_in) ) allocate ( thv_in(im,jm,lm_in) ) allocate ( dp_in(im,jm,lm_in) ) allocate ( pk_in(im,jm,lm_in) ) allocate ( ple_in(im,jm,lm_in+1) ) allocate ( ps_in(im,jm) ) allocate ( ak_in(lm_in+1) ) allocate ( bk_in(lm_in+1) ) read (10) ak_in read (10) bk_in do L=1,lm_in read(10) (( u_in(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm_in read(10) (( v_in(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm_in read(10) ((thv_in(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm_in+1 read(10) ((ple_in(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm_in read(10) (( pk_in(i,j,L),i=1,im),j=1,jm) enddo close (10) ps_in(:,:) = ple_in(:,:,lm_in+1) do L=lm_in,1,-1 dp_in(:,:,L) = ple_in(:,:,L+1) - ple_in(:,:,L) enddo C ********************************************************************** C **** Read Moist Internal Restart **** C ********************************************************************** allocate ( dum(im,jm) ) open (10,file=trim(moistrst),form='unformatted',access='sequential') nt = 0 rc = 0 dowhile (rc.eq.0) read (10,iostat=rc) dum if( rc.eq.0 ) then nt = nt + 1 allocate( q_out(im,jm,lm_in,nt) ) q_out(:,:,1,nt) = dum do L=2,lm_in read (10,iostat=rc) dum q_out(:,:,L,nt) = dum enddo if( nt.eq.1) then allocate( q_in (im,jm,lm_in,nt) ) q_in (:,:,:,nt) = q_out(:,:,:,nt) else q_out(:,:,:,1:nt-1) = q_in(:,:,:,1:nt-1) deallocate( q_in ) allocate( q_in (im,jm,lm_in,nt) ) q_in = q_out endif print *, 'Reading Moist Restart for Field # ',nt deallocate( q_out ) endif enddo close (10) ! ********************************************************************** ! **** Read Topography Datasets **** ! ********************************************************************** allocate ( phis(im,jm) ) print *, 'Reading Topography Dataset: ',trim(topo) open (10,file=trim(topo),form='unformatted',access='sequential') read (10) dum close(10) grav = getcon('GRAVITY') phis = dum*grav C ********************************************************************** C **** Remap State **** C ********************************************************************** allocate ( u_out(im,jm,lm_out) ) allocate ( v_out(im,jm,lm_out) ) allocate ( thv_out(im,jm,lm_out) ) allocate ( dp_out(im,jm,lm_out) ) allocate ( pk_out(im,jm,lm_out) ) allocate ( pke_out(im,jm,lm_out+1) ) allocate ( ple_out(im,jm,lm_out+1) ) allocate ( q_out(im,jm,lm_out,nt) ) allocate ( ps_out(im,jm) ) allocate ( ak_out(lm_out+1) ) allocate ( bk_out(lm_out+1) ) print *, 'Calling REMAP ...' call remap ( ps_out,dp_out,u_out,v_out,thv_out,q_out,phis,lm_out, . ps_in ,dp_in ,u_in ,v_in ,thv_in ,q_in ,phis,lm_in , . im,jm,nt,ak_out,bk_out ) print *, ' REMAP Finished' kappa = getcon('KAPPA') ple_out(:,:,lm_out+1) = ps_out(:,:) do L=lm_out,1,-1 ple_out(:,:,L) = ple_out(:,:,L+1)-dp_out(:,:,L) enddo c Ensure top edge = ptop c ---------------------- dp_out(:,:,1) = ple_out(:,:,2)-ak_out(1) ple_out(:,:,1) = ple_out(:,:,2)-dp_out(:,:,1) pke_out(:,:,:) = ple_out(:,:,:)**kappa do L=1,lm_out pk_out(:,:,L) = ( pke_out(:,:,L+1)-pke_out(:,:,L) ) . / ( kappa*log(ple_out(:,:,L+1)/ple_out(:,:,L)) ) enddo C ********************************************************************** C **** Write dycore internal Restart **** C ********************************************************************** write(date,101) nymd write(hour,102) nhms/10000 write(cim ,103) im write(cjm ,103) jm write(clm ,102) lm_out 101 format(i8.8) 102 format(i2.2) 103 format(i3.3) dynrst = trim(dynrst) // '.r' // cim // 'x' // cjm // 'x' // clm // . '.e' // date // '_' // hour // 'z' moistrst = trim(moistrst) // '.r' // cim // 'x' // cjm // 'x' // clm // . '.e' // date // '_' // hour // 'z' print * print *, ' Creating: ',trim(dynrst) print *, ' Creating: ',trim(moistrst) print *, ' output resolution: ',im,jm,lm_out print *, ' date: ',nymd,nhms print * open (10,file=trim(dynrst),form='unformatted',access='sequential') headr2(3) = lm_out write(10) headr1 write(10) headr2 write(10) ak_out write(10) bk_out do L=1,lm_out write(10) (( u_out(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm_out write(10) (( v_out(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm_out write(10) ((thv_out(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm_out+1 write(10) ((ple_out(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm_out write(10) (( pk_out(i,j,L),i=1,im),j=1,jm) enddo close (10) C ********************************************************************** C **** Write moist internal Restart **** C ********************************************************************** open (10,file=trim(moistrst),form='unformatted',access='sequential') do n=1,nt do L=1,lm_out dum(:,:) = q_out(:,:,L,n) write(10) dum enddo enddo close (10) stop end subroutine remap ( ps1,dp1,u1,v1,thv1,q1,phis1,lm1, . ps2,dp2,u2,v2,thv2,q2,phis2,lm2,im,jm,nt,ak,bk ) c ******************************************************************************* c ***** ***** c ***** Program to remap Target analysis variables (ps2,dp2,u2,v2,t2,q2) ***** c ***** onto Model grid variables (ps1,dp1,u1,v1,thv1,q1). Program ***** c ***** allows for differenct topographies between Analysis and Model. ***** c ***** ***** c ******************************************************************************* implicit none integer im,jm,lm1,lm2,nt c Model variables c --------------- real dp1(im,jm,lm1) real u1(im,jm,lm1) real v1(im,jm,lm1) real thv1(im,jm,lm1) real q1(im,jm,lm1,nt) real ps1(im,jm) real phis1(im,jm) c Target analysis variables c ------------------------- real dp2(im,jm,lm2) real u2(im,jm,lm2) real v2(im,jm,lm2) real t2(im,jm,lm2) real thv2(im,jm,lm2) real q2(im,jm,lm2,nt) real ps2(im,jm) real phis2(im,jm) real ak(lm2+1) real bk(lm2+1) c Local variables c --------------- real pz(im,jm) real pe0(im,jm,lm1+1) real pe1(im,jm,lm1+1) real pe2(im,jm,lm2+1) real pk (im,jm,lm2 ) real pke0(im,jm,lm1+1) real pke1(im,jm,lm1+1) real pke2(im,jm,lm2+1) real phi2(im,jm,lm2+1) real ptop1(im,jm) real getcon,kappa,cp,pl,dum,dum1,dum2 real rgas,pref,tref,pkref,tstar,eps,rvap,grav integer i,j,L kappa = getcon('KAPPA') rgas = getcon('RGAS') rvap = getcon('RVAP') grav = getcon('GRAVITY') cp = getcon('CP') eps = rvap/rgas-1.0 c Compute edge-level pressures c ---------------------------- pe1(:,:,lm1+1) = ps1(:,:) do L=lm1,1,-1 pe1(:,:,L) = pe1(:,:,L+1)-dp1(:,:,L) enddo ptop1(:,:) = pe1(:,:,1) c Construct target analysis pressure variables c -------------------------------------------- pe2(:,:,lm2+1) = ps2(:,:) do L=lm2,1,-1 pe2(:,:,L) = pe2(:,:,L+1) - dp2(:,:,L) enddo do j=1,jm do i=1,im if( pe2(i,j,1).eq.0.0 ) pe2(i,j,1) = 0.01*pe2(i,j,2) ! Prohibit ptop = 0 enddo enddo pke2(:,:,:) = pe2(:,:,:)**kappa c Construct target virtual potential temperature c ---------------------------------------------- do L=1,lm2 pk(:,:,L) = ( pke2(:,:,L+1)-pke2(:,:,L) )/( kappa*log(pe2(:,:,L+1)/pe2(:,:,L)) ) enddo c Construct target analysis heights c --------------------------------- phi2(:,:,lm2+1) = phis2(:,:) do L=lm2,1,-1 phi2(:,:,L) = phi2(:,:,L+1) + cp*thv2(:,:,L)*( pke2(:,:,L+1)-pke2(:,:,L) ) enddo c Compute new surface pressure consistent with target surface pressure and topography c ----------------------------------------------------------------------------------- do j=1,jm do i=1,im L = lm2 do while ( phi2(i,j,L).lt.phis1(i,j) ) L = L-1 enddo ps1(i,j) = pe2(i,j,L+1)*( 1+(phi2(i,j,L+1)-phis1(i,j))/(cp*thv2(i,j,L)*pke2(i,j,L+1)) )**(1.0/kappa) enddo enddo c Construct model pressure variables using new surface pressure c ------------------------------------------------------------- print * call set_eta ( lm1,dum,ak,bk ) L = 1 dum1 = (ak(L)+bk(L)*100000.0)/100 write(6,1000) L,ak(L),bk(L),dum1 1000 format(1x,'L: ',i3,4x,'ak: ',f10.3,' bk: ',f10.8,4x,'pe: ',f8.3) do L=2,lm1+1 dum2 = (ak(L)+bk(L)*100000.0)/100 write(6,1001) L,ak(L),bk(L),dum2,dum2-dum1 1001 format(1x,'L: ',i3,4x,'ak: ',f10.3,' bk: ',f10.8,4x,'pe: ',f8.3,3x,'dp: ',f7.3) dum1 = dum2 enddo print * do L=1,lm1+1 do j=1,jm do i=1,im pe1 (i,j,L) = ak(L)+bk(L)*ps1(i,j) pke1(i,j,L) = pe1(i,j,L)**kappa enddo enddo enddo do L=1,lm1 do j=1,jm do i=1,im dp1(i,j,L) = pe1(i,j,L+1)-pe1(i,j,L) enddo enddo enddo c Map target analysis onto grid defined by new surface pressure c ------------------------------------------------------------- print *, 'Calling GMAP, LM_in : ',lm2 print *, ' LM_out: ',lm1 call gmap ( im,jm,nt, kappa, . lm2, pke2, pe2, u2, v2, thv2, q2, . lm1, pke1, pe1, u1, v1, thv1, q1 ) return end c****6***0*********0*********0*********0*********0*********0**********72 subroutine gmap(im, jm, nq, akap, & km, pk3d_m, pe3d_m, u_m, v_m, pt_m, q_m, & kn, pk3d_n, pe3d_n, u_n, v_n, pt_n, q_n ) c****6***0*********0*********0*********0*********0*********0**********72 implicit none integer im, jm integer km, kn, nq C Input: C original data km-level real u_m(im,jm,km) real v_m(im,jm,km) real pt_m(im,jm,km) real q_m(im,jm,km,nq) real pk3d_m(im,jm,km+1) real pe3d_m(im,jm,km+1) C Output: C New data (kn-level) real u_n(im,jm,kn) real v_n(im,jm,kn) real pt_n(im,jm,kn) real q_n(im,jm,kn,nq) real pk3d_n(im,jm,kn+1) real pe3d_n(im,jm,kn+1) c local (private) integer i, j, k, ic, n real pe1(im,km+1),pe2(im,kn+1) real pk1(im,km+1),pk2(im,kn+1) real dp1(im,km) ,dp2(im,kn) real u1(im,km) , u2(im,kn) real v1(im,km) , v2(im,kn) real t1(im,km) , t2(im,kn) real q1(im,km) , q2(im,kn) real ptop real akap real ple, pek, dak, bkh real undef real big parameter ( undef = 1.e15 ) parameter ( big = 1.e10 ) do 2000 j=1,jm c Copy original data to local 2D arrays. do k=1,km+1 do i=1,im pe1(i,k) = pe3d_m(i,j,k) pk1(i,k) = pk3d_m(i,j,k) enddo enddo do k=1,kn+1 do i=1,im pe2(i,k) = pe3d_n(i,j,k) pk2(i,k) = pk3d_n(i,j,k) enddo enddo do k=1,km do i=1,im dp1(i,k) = pk1(i,k+1)-pk1(i,k) u1(i,k) = u_m(i,j,k) v1(i,k) = v_m(i,j,k) t1(i,k) = pt_m(i,j,k) enddo enddo do k=1,kn do i=1,im dp2(i,k) = pk2(i,k+1)-pk2(i,k) enddo enddo c map pt c ------ call mappm ( km, pk1, dp1, t1, kn, pk2, dp2, t2, im, 1, 7 ) do k=1,km do i=1,im dp1(i,k) = pe1(i,k+1)-pe1(i,k) enddo enddo do k=1,kn do i=1,im dp2(i,k) = pe2(i,k+1)-pe2(i,k) enddo enddo c map u,v c ------- call mappm ( km, pe1, dp1, u1, kn, pe2, dp2, u2, im, -1, 7 ) call mappm ( km, pe1, dp1, v1, kn, pe2, dp2, v2, im, -1, 7 ) c map q c ------- do n=1,nq do k=1,km do i=1,im q1(i,k) = q_m(i,j,k,n) enddo enddo call mappm ( km, pe1, dp1, q1, kn, pe2, dp2, q2, im, 0, 7 ) do k=1,kn do i=1,im q_n(i,j,k,n) = q2(i,k) enddo enddo enddo do k=1,kn do i=1,im u_n(i,j,k) = u2(i,k) v_n(i,j,k) = v2(i,k) pt_n(i,j,k) = t2(i,k) enddo enddo 2000 continue return end C****6***0*********0*********0*********0*********0*********0**********72 subroutine mappm(km, pe1, dp1, q1, kn, pe2, dp2, q2, im, iv, kord) C****6***0*********0*********0*********0*********0*********0**********72 C IV = 0: constituents C IV = 1: potential temp C IV =-1: winds C C Mass flux preserving mapping: q1(im,km) -> q2(im,kn) C C pe1: pressure at layer edges (from model top to bottom surface) C in the original vertical coordinate C pe2: pressure at layer edges (from model top to bottom surface) C in the new vertical coordinate parameter (kmax = 200) parameter (R3 = 1./3., R23 = 2./3.) real dp1(im,km), dp2(im,kn), & q1(im,km), q2(im,kn), & pe1(im,km+1), pe2(im,kn+1) integer kord C local work arrays real a4(4,im,km) do k=1,km do i=1,im a4(1,i,k) = q1(i,k) enddo enddo call ppm2m(a4, dp1, im, km, iv, kord) C Lowest layer: constant distribution do i=1, im a4(2,i,km) = q1(i,km) a4(3,i,km) = q1(i,km) a4(4,i,km) = 0. enddo do 5555 i=1,im k0 = 1 do 555 k=1,kn if(pe2(i,k+1) .le. pe1(i,1)) then ! Entire grid above old ptop q2(i,k) = a4(2,i,1) elseif(pe2(i,k) .ge. pe1(i,km+1)) then ! Entire grid below old ps q2(i,k) = a4(3,i,km) elseif(pe2(i,k ) .lt. pe1(i,1) .and. & pe2(i,k+1) .gt. pe1(i,1)) then ! Part of the grid above ptop q2(i,k) = a4(1,i,1) else do 45 L=k0,km ! locate the top edge at pe2(i,k) if( pe2(i,k) .ge. pe1(i,L) .and. & pe2(i,k) .le. pe1(i,L+1) ) then k0 = L PL = (pe2(i,k)-pe1(i,L)) / dp1(i,L) if(pe2(i,k+1) .le. pe1(i,L+1)) then ! entire new grid is within the original grid PR = (pe2(i,k+1)-pe1(i,L)) / dp1(i,L) TT = R3*(PR*(PR+PL)+PL**2) q2(i,k) = a4(2,i,L) + 0.5*(a4(4,i,L)+a4(3,i,L) & - a4(2,i,L))*(PR+PL) - a4(4,i,L)*TT goto 555 else ! Fractional area... delp = pe1(i,L+1) - pe2(i,k) TT = R3*(1.+PL*(1.+PL)) qsum = delp*(a4(2,i,L)+0.5*(a4(4,i,L)+ & a4(3,i,L)-a4(2,i,L))*(1.+PL)-a4(4,i,L)*TT) dpsum = delp k1 = L + 1 goto 111 endif endif 45 continue 111 continue do 55 L=k1,km if( pe2(i,k+1) .gt. pe1(i,L+1) ) then ! Whole layer.. qsum = qsum + dp1(i,L)*q1(i,L) dpsum = dpsum + dp1(i,L) else delp = pe2(i,k+1)-pe1(i,L) esl = delp / dp1(i,L) qsum = qsum + delp * (a4(2,i,L)+0.5*esl* & (a4(3,i,L)-a4(2,i,L)+a4(4,i,L)*(1.-r23*esl)) ) dpsum = dpsum + delp k0 = L goto 123 endif 55 continue delp = pe2(i,k+1) - pe1(i,km+1) if(delp .gt. 0.) then ! Extended below old ps qsum = qsum + delp * a4(3,i,km) dpsum = dpsum + delp endif 123 q2(i,k) = qsum / dpsum endif 555 continue 5555 continue return end c****6***0*********0*********0*********0*********0*********0**********72 subroutine ppm2m(a4,delp,im,km,iv,kord) c****6***0*********0*********0*********0*********0*********0**********72 c iv = 0: positive definite scalars c iv = 1: others c iv =-1: winds implicit none integer im, km, lmt, iv integer kord integer i, k, km1 real a4(4,im,km), delp(im,km) c local arrays. real dc(im,km),delq(im,km) real h2(im,km) real a1, a2, a3, b2, c1, c2, c3, d1, d2, f1, f2, f3, f4 real s1, s2, s3, s4, ss3, s32, s34, s42, sc real qmax, qmin, cmax, cmin real dm, qm, dq, tmp C Local scalars: real qmp real lac km1 = km - 1 do 500 k=2,km do 500 i=1,im delq(i,k-1) = a4(1,i,k) - a4(1,i,k-1) 500 a4(4,i,k ) = delp(i,k-1) + delp(i,k) do 1220 k=2,km1 do 1220 i=1,im c1 = (delp(i,k-1)+0.5*delp(i,k))/a4(4,i,k+1) c2 = (delp(i,k+1)+0.5*delp(i,k))/a4(4,i,k) tmp = delp(i,k)*(c1*delq(i,k) + c2*delq(i,k-1)) / & (a4(4,i,k)+delp(i,k+1)) qmax = max(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1)) - a4(1,i,k) qmin = a4(1,i,k) - min(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1)) dc(i,k) = sign(min(abs(tmp),qmax,qmin), tmp) 1220 continue c****6***0*********0*********0*********0*********0*********0**********72 c 4th order interpolation of the provisional cell edge value c****6***0*********0*********0*********0*********0*********0**********72 do 12 k=3,km1 do 12 i=1,im c1 = delq(i,k-1)*delp(i,k-1) / a4(4,i,k) a1 = a4(4,i,k-1) / (a4(4,i,k) + delp(i,k-1)) a2 = a4(4,i,k+1) / (a4(4,i,k) + delp(i,k)) a4(2,i,k) = a4(1,i,k-1) + c1 + 2./(a4(4,i,k-1)+a4(4,i,k+1)) * & ( delp(i,k)*(c1*(a1 - a2)+a2*dc(i,k-1)) - & delp(i,k-1)*a1*dc(i,k ) ) 12 continue C Area preserving cubic with 2nd deriv. = 0 at the boundaries C Top do i=1,im d1 = delp(i,1) d2 = delp(i,2) qm = (d2*a4(1,i,1)+d1*a4(1,i,2)) / (d1+d2) dq = 2.*(a4(1,i,2)-a4(1,i,1)) / (d1+d2) c1 = 4.*(a4(2,i,3)-qm-d2*dq) / ( d2*(2.*d2*d2+d1*(d2+3.*d1)) ) c3 = dq - 0.5*c1*(d2*(5.*d1+d2)-3.*d1**2) a4(2,i,2) = qm - 0.25*c1*d1*d2*(d2+3.*d1) a4(2,i,1) = d1*(2.*c1*d1**2-c3) + a4(2,i,2) dc(i,1) = a4(1,i,1) - a4(2,i,1) C No over- and undershoot condition cmax = max(a4(1,i,1), a4(1,i,2)) cmin = min(a4(1,i,1), a4(1,i,2)) a4(2,i,2) = max(cmin,a4(2,i,2)) a4(2,i,2) = min(cmax,a4(2,i,2)) enddo if(iv == 0) then do i=1,im a4(2,i,1) = max(0.,a4(2,i,1)) a4(2,i,2) = max(0.,a4(2,i,2)) enddo elseif(iv == -1) then do i=1,im if( a4(2,i,1)*a4(1,i,1) <= 0. ) a4(2,i,1) = 0. enddo endif c****6***0*********0*********0*********0*********0*********0**********72 c Bottom c Area preserving cubic with 2nd deriv. = 0 at the surface do 15 i=1,im d1 = delp(i,km) d2 = delp(i,km1) qm = (d2*a4(1,i,km)+d1*a4(1,i,km1)) / (d1+d2) dq = 2.*(a4(1,i,km1)-a4(1,i,km)) / (d1+d2) c1 = (a4(2,i,km1)-qm-d2*dq) / (d2*(2.*d2*d2+d1*(d2+3.*d1))) c3 = dq - 2.0*c1*(d2*(5.*d1+d2)-3.*d1**2) a4(2,i,km) = qm - c1*d1*d2*(d2+3.*d1) a4(3,i,km) = d1*(8.*c1*d1**2-c3) + a4(2,i,km) dc(i,km) = a4(3,i,km) - a4(1,i,km) c****6***0*********0*********0*********0*********0*********0**********72 c No over- and undershoot condition cmax = max(a4(1,i,km), a4(1,i,km1)) cmin = min(a4(1,i,km), a4(1,i,km1)) a4(2,i,km) = max(cmin,a4(2,i,km)) a4(2,i,km) = min(cmax,a4(2,i,km)) c****6***0*********0*********0*********0*********0*********0**********72 15 continue if(iv .eq. 0) then do i=1,im a4(2,i,km) = max(0.,a4(2,i,km)) a4(3,i,km) = max(0.,a4(3,i,km)) enddo endif do 20 k=1,km1 do 20 i=1,im a4(3,i,k) = a4(2,i,k+1) 20 continue c c f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 ) c c Top 2 and bottom 2 layers always use monotonic mapping do k=1,2 do i=1,im a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) enddo call kmppm(dc(1,k),a4(1,1,k),im, 0) enddo if(kord == 7) then c****6***0*********0*********0*********0*********0*********0**********72 C Huynh's 2nd constraint c****6***0*********0*********0*********0*********0*********0**********72 do k=2, km1 do i=1,im h2(i,k) = delq(i,k) - delq(i,k-1) enddo enddo do 4000 k=3, km-2 do 3000 i=1, im C Right edges qmp = a4(1,i,k) + 2.0*delq(i,k-1) lac = a4(1,i,k) + 1.5*h2(i,k-1) + 0.5*delq(i,k-1) qmin = min(a4(1,i,k), qmp, lac) qmax = max(a4(1,i,k), qmp, lac) a4(3,i,k) = min(max(a4(3,i,k), qmin), qmax) C Left edges qmp = a4(1,i,k) - 2.0*delq(i,k) lac = a4(1,i,k) + 1.5*h2(i,k+1) - 0.5*delq(i,k) qmin = min(a4(1,i,k), qmp, lac) qmax = max(a4(1,i,k), qmp, lac) a4(2,i,k) = min(max(a4(2,i,k), qmin), qmax) C Recompute A6 a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) 3000 continue ! Additional constraint to prevent negatives if (iv == 0) then call kmppm(dc(1,k),a4(1,1,k),im, 2) endif 4000 continue else lmt = kord - 3 lmt = max(0, lmt) if (iv .eq. 0) lmt = min(2, lmt) do k=3, km-2 do i=1,im a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) enddo call kmppm(dc(1,k),a4(1,1,k),im, lmt) enddo endif do 5000 k=km1,km do i=1,im a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) enddo call kmppm(dc(1,k),a4(1,1,k),im, 0) 5000 continue return end c****6***0*********0*********0*********0*********0*********0**********72 subroutine kmppm(dm, a4, km, lmt) c****6***0*********0*********0*********0*********0*********0**********72 implicit none real r12 parameter (r12 = 1./12.) integer km, lmt integer i real a4(4,km),dm(km) real da1, da2, a6da real fmin real qmp if (lmt .eq. 3) return ! Full constraint if(lmt .eq. 0) then do 100 i=1,km if(dm(i) .eq. 0.) then a4(2,i) = a4(1,i) a4(3,i) = a4(1,i) a4(4,i) = 0. else da1 = a4(3,i) - a4(2,i) da2 = da1**2 a6da = a4(4,i)*da1 if(a6da .lt. -da2) then a4(4,i) = 3.*(a4(2,i)-a4(1,i)) a4(3,i) = a4(2,i) - a4(4,i) elseif(a6da .gt. da2) then a4(4,i) = 3.*(a4(3,i)-a4(1,i)) a4(2,i) = a4(3,i) - a4(4,i) endif endif 100 continue elseif (lmt .eq. 2) then c Positive definite c Positive definite constraint do 250 i=1,km if(abs(a4(3,i)-a4(2,i)) .ge. -a4(4,i)) go to 250 fmin = a4(1,i)+0.25*(a4(3,i)-a4(2,i))**2/a4(4,i)+a4(4,i)*r12 if(fmin.ge.0.) go to 250 if(a4(1,i).lt.a4(3,i) .and. a4(1,i).lt.a4(2,i)) then a4(3,i) = a4(1,i) a4(2,i) = a4(1,i) a4(4,i) = 0. elseif(a4(3,i) .gt. a4(2,i)) then a4(4,i) = 3.*(a4(2,i)-a4(1,i)) a4(3,i) = a4(2,i) - a4(4,i) else a4(4,i) = 3.*(a4(3,i)-a4(1,i)) a4(2,i) = a4(3,i) - a4(4,i) endif 250 continue elseif (lmt == 1) then ! Improved full monotonicity constraint (Lin) ! Note: no need to provide first guess of A6 <-- a4(4,i) do i=1, km qmp = 2.*dm(i) a4(2,i) = a4(1,i)-sign(min(abs(qmp),abs(a4(2,i)-a4(1,i))), qmp) a4(3,i) = a4(1,i)+sign(min(abs(qmp),abs(a4(3,i)-a4(1,i))), qmp) a4(4,i) = 3.*( 2.*a4(1,i) - (a4(2,i)+a4(3,i)) ) enddo endif return end subroutine minmax ( q,name,im,jm,lm ) character*4 name real q(im,jm,lm) qmax = q(1,1,1) qmin = q(1,1,1) do l=1,lm do j=1,jm do i=1,im qmax = max( qmax,q(i,j,L) ) qmin = min( qmin,q(i,j,L) ) enddo enddo enddo print *, name,' max: ',qmax,' min: ',qmin return end subroutine usage() print *, "Usage: " print * print *, " rs_vinterp.x -dynrst dynrst_fname Default: fvcore_internal_restart" print *, " -moistrst moistrst_fname Default: moist_internal_restart" print *, " -topo topo_fname" print *, " -lm LM" print * print *, "where:" print *, "-----" print *, " -dynrst dynrst_fname: Filename of dynamics internal restart" print *, " -moistrst moistrst_fname: Filename of moist internal restart" print *, " -topo topo_fname: Filename of topography" print *, " -lm LM: Output Veritcal Resolution" print * print *, "creates updated restarts at new LM resolution" print *, "---------------------------------------------" print * call exit(7) end FUNCTION GETCON (NAME) C*********************************************************************** C C GENERIC FUNCTION GETCON IS A REPOSITORY OF GLOBAL VARIABLES, I.E. C A MEMORY FOR SCALAR VALUES NEEDED THROUGHOUT A LARGE PROGRAM. C C THIS SPECIFIC FUNCTION, GETCON, REMEMBERS FLOATING POINT VALUES. C THE FUNCTION IS CALLED WITH A CHARACTER NAME TO INTERROGATE A VALUE. C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** PARAMETER (MAXCON=40) CHARACTER*(*) NAME CHARACTER*16 ANAME(MAXCON) real ACON (MAXCON) C COMPUTATIONAL CONSTANTS C ----------------------- PARAMETER ( VECMAX = 65535.5 ) PARAMETER ( CALTOJ = 4184. ) PARAMETER ( UNDEF = 1.E15 ) C ASTRONOMICAL CONSTANTS C ---------------------- PARAMETER ( OB = 23.45 ) PARAMETER ( PER = 102. ) PARAMETER ( ECC = .0167 ) PARAMETER ( AE = 6371E3 ) PARAMETER ( EQNX = 80.5 ) PARAMETER ( SOLS = 176.5 ) PARAMETER ( S0 = 1365.0 ) C TERRESTRIAL CONSTANTS C --------------------- c PARAMETER ( GRAV = 9.81 ) ! GEOS PARAMETER ( GRAV = 9.80616 ) ! fv PARAMETER ( SRFPRS = 984.7 ) PARAMETER ( PIMEAN = 984.7 ) PARAMETER ( PSTD = 1000.0 ) PARAMETER ( TSTD = 280.0 ) PARAMETER ( SDAY = 86400.0 ) PARAMETER ( SSALB = 0.99 ) PARAMETER ( CO2 = 355.0 ) PARAMETER ( CFC11 = 0.3 ) PARAMETER ( CFC12 = 0.5 ) PARAMETER ( CFC22 = 0.2 ) C THERMODYNAMIC CONSTANTS C ----------------------- c PARAMETER ( CPD = 1004.16 ) ! GEOS PARAMETER ( CPD = 1004.64 ) ! fv PARAMETER ( CPV = 1869.46 ) PARAMETER ( ALHL = 2.499E6 ) PARAMETER ( ALHS = 2.845E6 ) PARAMETER ( STFBOL = 5.67E-8 ) PARAMETER ( AIRMW = 28.97 ) PARAMETER ( H2OMW = 18.01 ) PARAMETER ( RUNIV = 8314.3 ) c PARAMETER ( RGAS = RUNIV/AIRMW) ! GEOS c PARAMETER ( RVAP = RUNIV/H2OMW) ! GEOS PARAMETER ( RGAS = 287.04 ) ! fv PARAMETER ( RVAP = 461.00 ) ! fv PARAMETER ( RKAP = RGAS/CPD ) PARAMETER ( HEATW = 597.2 ) PARAMETER ( HEATI = 680.0 ) PARAMETER ( TICE = 273.16 ) C TURBULENCE CONSTANTS C -------------------- PARAMETER ( VKRM = 0.4 ) C MOISTURE CONSTANTS C ------------------ PARAMETER ( EPS = 0.622 ) PARAMETER ( VIRTCON= 0.609 ) PARAMETER ( EPSFAC = EPS*HEATW/RGAS*CALTOJ ) DATA ANAME(1 ),ACON(1 ) / 'CP ', CPD / DATA ANAME(2 ),ACON(2 ) / 'RGAS ', RGAS / DATA ANAME(3 ),ACON(3 ) / 'KAPPA ', RKAP / DATA ANAME(4 ),ACON(4 ) / 'LATENT HEAT COND', ALHL / DATA ANAME(5 ),ACON(5 ) / 'GRAVITY ', GRAV / DATA ANAME(6 ),ACON(6 ) / 'STEFAN-BOLTZMAN ', STFBOL / DATA ANAME(7 ),ACON(7 ) / 'VON KARMAN ', VKRM / DATA ANAME(8 ),ACON(8 ) / 'EARTH RADIUS ', AE / DATA ANAME(9 ),ACON(9 ) / 'OBLIQUITY ', OB / DATA ANAME(10),ACON(10) / 'ECCENTRICITY ', ECC / DATA ANAME(11),ACON(11) / 'PERIHELION ', PER / DATA ANAME(12),ACON(12) / 'VERNAL EQUINOX ', EQNX / DATA ANAME(13),ACON(13) / 'SUMMER SOLSTICE ', SOLS / DATA ANAME(14),ACON(14) / 'MAX VECT LENGTH ', VECMAX / DATA ANAME(15),ACON(15) / 'MOL WT H2O ', H2OMW / DATA ANAME(16),ACON(16) / 'MOL WT AIR ', AIRMW / DATA ANAME(17),ACON(17) / 'CPV ', CPV / DATA ANAME(18),ACON(18) / 'CPD ', CPD / DATA ANAME(19),ACON(19) / 'UNIV GAS CONST ', RUNIV / DATA ANAME(20),ACON(20) / 'LATENT HEAT SBLM', ALHS / DATA ANAME(21),ACON(21) / 'FREEZING-POINT ', TICE / DATA ANAME(23),ACON(23) / 'CALTOJ ', CALTOJ / DATA ANAME(24),ACON(24) / 'EPS ', EPS / DATA ANAME(25),ACON(25) / 'HEATW ', HEATW / DATA ANAME(26),ACON(26) / 'EPSFAC ', EPSFAC / DATA ANAME(27),ACON(27) / 'VIRTCON ', VIRTCON/ DATA ANAME(28),ACON(28) / 'PIMEAN ', PIMEAN / DATA ANAME(29),ACON(29) / 'SDAY ', SDAY / DATA ANAME(30),ACON(30) / 'HEATI ', HEATI / DATA ANAME(31),ACON(31) / 'S0 ', S0 / DATA ANAME(32),ACON(32) / 'PSTD ', PSTD / DATA ANAME(33),ACON(33) / 'TSTD ', TSTD / DATA ANAME(34),ACON(34) / 'SSALB ', SSALB / DATA ANAME(35),ACON(35) / 'UNDEF ', UNDEF / DATA ANAME(36),ACON(36) / 'CO2 ', CO2 / DATA ANAME(37),ACON(37) / 'RVAP ', RVAP / DATA ANAME(38),ACON(38) / 'CFC11 ', CFC11 / DATA ANAME(39),ACON(39) / 'CFC12 ', CFC12 / DATA ANAME(40),ACON(40) / 'CFC22 ', CFC22 / DO 10 I=1,MAXCON IF(NAME.EQ.ANAME(I)) THEN GETCON = ACON(I) RETURN ENDIF 10 CONTINUE 900 PRINT *,' CANNOT FIND FLOATING POINT CONSTANT - ',NAME PRINT *,' GETCON - CANNOT FIND CONSTANT REQUESTED' END