<p><b>fanglin.yang@noaa.gov</b> 2012-05-09 19:07:10 -0600 (Wed, 09 May 2012)</p><p>9may2012 work version<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics_gfs/src/core_atmos_physics_gfs/physics_gfs/dotstep.f
===================================================================
--- branches/atmos_physics_gfs/src/core_atmos_physics_gfs/physics_gfs/dotstep.f 2012-05-09 22:12:47 UTC (rev 1884)
+++ branches/atmos_physics_gfs/src/core_atmos_physics_gfs/physics_gfs/dotstep.f 2012-05-10 01:07:10 UTC (rev 1885)
@@ -58,11 +58,11 @@
real(kind=kind_phys) :: sfc_mpas(nsfc_mpas,ncell_mpas)
real(kind=kind_phys) :: air_mpas(nair_mpas,ncell_mpas,levs_mpas)
-! --prsi : model interface level pressure in Pasca
-! --prsl : model integer layer pressure in Pasca
+! --prsi : model interface level pressure in centibar
+! --prsl : model integer layer pressure in centibar
! --gu : model layer zonal wind m/s
! --gv : model layer meridional wind m/s
-! --vvel : model layer vertical velocity in Pascal/sec
+! --vvel : model layer vertical velocity in centibar/sec
! --gt : model layer temperature in K
! --gr : model layer specific humidity in gm/gm
! --gtrace: model layer tracer (ozne and cloud water) mass mixing ratio
@@ -72,23 +72,17 @@
!-----mpas related fileds--------
!****************************************************************************
- INTEGER :: INDLSEV,JBASEV, INDLSOD,JBASOD
- include 'function2'
-!
-!
-!------------------------------------------------------------
+
!--at initial time to allocate GFS-related model arrays,
!--to read in GFS namelist file and set up running parameters,
!--and to read in most static boundary conditions.
-
-
if(ifirst > 0) then
!--GFS initial condition date, 1-hour,2-month,3-day,4-year
do j=1,4
idate(j)=idate_mpas(j)
enddo
- lats_node_r=1 !for MPAS use 1-D block for each node
+ lats_node_r=1 !for MPAS use 1-D block for each task
call GFS_Initialize(node0_mpas,fhour_mpas,levs_mpas,ncell_mpas,
& xlon_mpas,xlat_mpas,lats_node_r,dt_mpas,nlunit_mpas,
@@ -179,7 +173,7 @@
if(.not. adiab) then
call gloopr
!---input
- & (lonsperlar,global_lats_r,phour,xlon,xlat,kdt,
+ & (phour,kdt,lonsperlar,global_lats_r,xlon,xlat,
& sfc_fld%slmsk,sfc_fld%sheleg,
& sfc_fld%zorl, sfc_fld%tsea,
& sfc_fld%alvsf, sfc_fld%alnsf, sfc_fld%alvwf,
@@ -200,18 +194,15 @@
endif !sswr .or. lslwr
-!set to zero for every timestep
-
if(.not. adiab) then
call gloopb
- & (lonsperlar,
- & deltim,phour,sfc_fld, flx_fld, nst_fld, SFALB,
- & xlon,
+ & (phour,kdt,deltim,lonsperlar,global_lats_r,
+ & lsout,fscav,xlon,xlat,
+ & sfc_fld, flx_fld, nst_fld, sfalb,
& swh,hlw,hprime,slag,sdec,cdec,
- & ozplin,jindx1,jindx2,ddy,pdryini,
- & phy_f3d, phy_f2d, xlat,kdt,
- & batah,lsout,fscav,ngptc)
- endif ! not.adiab
+ & ozplin,jindx1,jindx2,ddy,
+ & phy_f3d, phy_f2d)
+ endif
!
IF (mod(kdt,nszer) == 0 .and. lsout) THEN
Added: branches/atmos_physics_gfs/src/core_atmos_physics_gfs/physics_gfs/gloopb.f
===================================================================
--- branches/atmos_physics_gfs/src/core_atmos_physics_gfs/physics_gfs/gloopb.f (rev 0)
+++ branches/atmos_physics_gfs/src/core_atmos_physics_gfs/physics_gfs/gloopb.f 2012-05-10 01:07:10 UTC (rev 1885)
@@ -0,0 +1,726 @@
+!--modeified for MPAS, Fanglin Yang, May 2012
+!
+ subroutine gloopb
+ & (phour,kdt,tstep,lonsperlar, global_lats_r,
+ & lsout,fscav,xlon,xlat,
+ & sfc_fld, flx_fld, nst_fld, sfalb,
+ & swh,hlw,hprime,slag,sdec,cdec,
+ & ozplin,jindx1,jindx2,ddy,
+ & phy_f3d, phy_f2d)
+!!
+#include "f_hpm.h"
+!!
+ use machine , only : kind_evod,kind_phys,kind_rad
+ use resol_def , only : jcap,latr,levs,lonr,
+ & lsoil,ncld,nmtvr,nrcm,ntcw,ntoz,
+ & ntrac,num_p2d,num_p3d,
+ & thermodyn_id,sfcpress_id,nfxr
+
+ use layout1 , only : ipt_lats_node_r,
+ & lat1s_r,lats_dim_r,
+ & lats_node_a,lats_node_r,
+ & me,nodes
+ use gg_def , only : coslat_r,rcs2_r,sinlat_r,wgt_r
+ use date_def , only : fhour,idate
+ use namelist_def , only : crtrh,fhswr,flgmin,
+ & gen_coord_hybrid,gg_tracers,
+ & hybrid,ldiag3d,lscca,lsfwd,
+ & lsm,lssav,lsswr,ncw,ngptc,
+ & old_monin,pre_rad,random_clds,
+ & ras,shuff_lats_r,
+ & sashal,ctei_rm,mom4ice,newsas,
+ & ccwf,cnvgwd,lggfs3d,trans_trac,
+ & mstrat,cal_pre,nst_fcst,
+ & dlqf,moist_adj,cdmbgwd,
+ & bkgd_vdif_m, bkgd_vdif_h,
+ & bkgd_vdif_s,shal_cnv,
+ & psautco, prautco, evpco, wminco
+ use coordinate_def , only : vertcoord_id
+ use module_ras , only : ras_init
+ use physcons , only : grav => con_g,
+ & rerth => con_rerth, ! hmhj
+ & fv => con_fvirt, ! mjr
+ & rvrdm1 => con_FVirt,
+ & rd => con_rd
+ use ozne_def , only : latsozp,levozp,
+ & pl_coeff,pl_pres,timeoz
+
+ use Sfc_Flx_ESMFMod
+ use Nst_Var_ESMFMod
+ use mersenne_twister
+ use d3d_def
+ use tracer_const
+!
+ include 'mpif.h'
+ implicit none
+!
+ TYPE(Sfc_Var_Data) :: sfc_fld
+ TYPE(Flx_Var_Data) :: flx_fld
+ TYPE(Nst_Var_Data) :: nst_fld
+!
+ real(kind=kind_phys), PARAMETER :: RLAPSE=0.65E-2
+ real(kind=kind_evod), parameter :: cons_0=0.0, cons_24=24.0
+ &, cons_99=99.0, cons_1p0d9=1.0E9
+
+
+!$$$ integer n1rac, n2rac,nlons_v(ngptc)
+!$$$ parameter (n1rac=ntrac-ntshft-1, n2rac=n1rac+1)
+!
+! integer id,njeff,istrt,lon,kdt
+ integer id,njeff, lon,kdt
+!!
+ real(kind=kind_phys) prsl(ngptc,levs)
+ real(kind=kind_phys) prslk(ngptc,levs),dpshc(ngptc)
+ real(kind=kind_phys) prsi(ngptc,levs+1),phii(ngptc,levs+1)
+ real(kind=kind_phys) prsik(ngptc,levs+1),phil(ngptc,levs)
+!!
+ real (kind=kind_phys) gu(ngptc,levs), gv1(ngptc,levs)
+ real (kind=kind_phys) ugrd(ngptc,levs),vgrd(ngptc,levs)
+ real (kind=kind_phys) gphi(ngptc), glam(ngptc)
+ real (kind=kind_phys) gq(ngptc), gt(ngptc,levs), pgr(ngptc)
+ real (kind=kind_phys) gr(ngptc,levs,ntrac)
+ real (kind=kind_phys) gd(ngptc,levs)
+ real (kind=kind_phys) adt(ngptc,levs), adr(ngptc,levs,ntrac)
+ real (kind=kind_phys) adu(ngptc,levs), adv(ngptc,levs)
+ real (kind=kind_phys) gtv(ngptc,levs) ! hmhj
+ real (kind=kind_phys) gtvx(ngptc,levs), gtvy(ngptc,levs) ! hmhj
+ real (kind=kind_phys) sumq(ngptc,levs), xcp(ngptc,levs)
+!
+ real (kind=kind_phys) dt3dt_v(ngptc,levs,6),
+ & dq3dt_v(ngptc,levs,5+pl_coeff),
+ & du3dt_v(ngptc,levs,4),
+ & dv3dt_v(ngptc,levs,4)
+ &, upd_mf_v(ngptc,levs), dwn_mf_v(ngptc,levs)
+ &, det_mf_v(ngptc,levs)
+ &, dkh_v(ngptc,LEVS), rnp_v(ngptc,levs)
+!!
+ real(kind=kind_evod) gq_save(lonr,lats_dim_r)
+!!
+ real (kind=kind_rad) slag,sdec,cdec,phour
+ real (kind=kind_rad) xlon(lonr,lats_node_r)
+ real (kind=kind_rad) xlat(lonr,lats_node_r)
+ real (kind=kind_rad) coszdg(lonr,lats_node_r),
+ & hprime(lonr,nmtvr,lats_node_r),
+! & fluxr(lonr,nfxr,lats_node_r),
+ & sfalb(lonr,lats_node_r)
+ real (kind=kind_rad) swh(lonr,levs,lats_node_r)
+ real (kind=kind_rad) hlw(lonr,levs,lats_node_r)
+!!
+ real (kind=kind_phys)
+ & phy_f3d(lonr,levs,num_p3d,lats_node_r),
+ & phy_f2d(lonr,num_p2d,lats_node_r), fscav(ntrac-ncld-1)
+!
+!
+ real (kind=kind_phys) exp,dtphys,dtp,dtf,sumed(2)
+ real (kind=kind_evod) tstep
+ real (kind=kind_phys) pdryini,sigshc,rk
+!!
+ integer global_lats_r(latr)
+ integer lonsperlar(latr)
+ integer dimg
+cc
+ integer i,ierr,iter,j,k,kap,kar,kat,kau,kav,ksq,jj,kk
+ integer kst,kdtphi,kdtlam ! hmhj
+ integer l,lan,lan0,lat,lmax,locl,ii,lonrbm
+! integer lon_dim,lons_lat,n,node
+ integer lons_lat,n,node
+ integer nsphys
+!
+ real(kind=kind_evod) pwatg(latr),pwatj(lats_node_r),
+ & pwatp,ptotg(latr),sumwa,sumto,
+ & ptotj(lats_node_r),pcorr,pdryg,
+ & solhr,clstp
+cc
+ integer ipt_ls ! hmhj
+ real(kind=kind_evod) reall ! hmhj
+
+ real(kind=kind_evod) typical_pgr
+c
+cc
+ integer indlsev,jbasev,n0
+ integer indlsod,jbasod
+cc
+ include 'function2'
+cc
+ real(kind=kind_evod) cons0,cons2 !constant
+cc
+ logical lsout
+ logical, parameter :: flipv = .true.
+cc
+! for nasa/nrl ozone production and distruction rates:(input through fixio)
+ real ozplin(latsozp,levozp,pl_coeff,timeoz)
+ integer jindx1(lats_node_r),jindx2(lats_node_r)!for ozone interpolaton
+ real ddy(lats_node_r) !for ozone interpolaton
+ real ozplout(levozp,lats_node_r,pl_coeff)
+!Moor real ozplout(lonr,levozp,pl_coeff,lats_node_r)
+!!
+ real(kind=kind_phys), allocatable :: acv(:,:),acvb(:,:),acvt(:,:)
+ save acv,acvb,acvt
+!!
+! integer, parameter :: maxran=5000
+! integer, parameter :: maxran=3000
+ integer, parameter :: maxran=6000, maxsub=6, maxrs=maxran/maxsub
+ type (random_stat) :: stat(maxrs)
+ real (kind=kind_phys), allocatable, save :: rannum_tank(:,:,:)
+ real (kind=kind_phys) :: rannum(lonr*latr)
+ integer iseed, nrc, seed0, kss, ksr, indxr(nrcm), iseedl
+ integer nf0,nf1,ind,nt,indod,indev
+ real(kind=kind_evod) fd2, wrk(1), wrk2(nrcm)
+
+ logical first,ladj
+ parameter (ladj=.true.)
+ data first/.true./
+! save krsize, first, nrnd,seed0
+ save first, seed0
+!!
+ integer nlons_v(ngptc)
+ real(kind=kind_phys) smc_v(ngptc,lsoil),stc_v(ngptc,lsoil)
+ &, slc_v(ngptc,lsoil)
+ &, swh_v(ngptc,levs), hlw_v(ngptc,levs)
+ &, vvel(ngptc,levs)
+ &, hprime_v(ngptc,nmtvr)
+ real(kind=kind_phys) phy_f3dv(ngptc,LEVS,num_p3d),
+ & phy_f2dv(ngptc,num_p2d)
+ &, rannum_v(ngptc,nrcm)
+ real(kind=kind_phys) sinlat_v(ngptc),coslat_v(ngptc)
+ &, ozplout_v(ngptc,levozp,pl_coeff)
+ real (kind=kind_rad) rqtk(ngptc), rcs2_lan, rcs_lan
+! real (kind=kind_rad) rcs_v(ngptc), rqtk(ngptc), rcs2_lan
+!
+!--------------------------------------------------------------------
+! print *,' in gloopb vertcoord_id =',vertcoord_id
+
+! real(kind=kind_evod) sinlat_v(lonr),coslat_v(lonr),rcs2_v(lonr)
+! real(kind=kind_phys) dpshc(lonr)
+ real (kind=kind_rad) work1, qmin, tem
+ parameter (qmin=1.0e-10)
+ integer ksd,ksplam,kspphi
+ integer ksu,ksv,ksz,item,jtem,ktem,ltem,mtem
+
+!
+! --- for debug test use
+ real (kind=kind_phys) :: temlon, temlat, alon, alat
+ integer :: ipt
+ logical :: lprnt
+!
+!
+!
+!----------------------
+ if (first) then
+!----------------------
+
+ allocate (acv(lonr,lats_node_r))
+ allocate (acvb(lonr,lats_node_r))
+ allocate (acvt(lonr,lats_node_r))
+!
+ seed0 = idate(1) + idate(2) + idate(3) + idate(4)
+ call random_setseed(seed0)
+ call random_number(wrk)
+ seed0 = seed0 + nint(wrk(1)*1000.0)
+!
+ if (.not. newsas) then ! random number needed for RAS and old SAS
+ if (random_clds) then ! create random number tank
+! -------------------------
+ if (.not. allocated(rannum_tank))
+ & allocate (rannum_tank(lonr,maxran,lats_node_r))
+! lonrb2 = lonr / 2
+ lonrbm = lonr / maxsub
+ if (me == 0) write(0,*)' maxran=',maxran,' maxrs=',maxrs,
+ & 'maxsub=',maxsub,' lonrbm=',lonrbm
+!$OMP parallel do private(nrc,iseedl,rannum,lat,i,j,k,ii,jj,kk)
+ do nrc=1,maxrs
+ iseedl = seed0 + nrc - 1
+ call random_setseed(iseedl,stat(nrc))
+ call random_number(rannum,stat(nrc))
+ do j=1,lats_node_r
+ lat = global_lats_r(ipt_lats_node_r-1+j)
+ jj = (lat-1)*lonr
+ do k=1,maxsub
+ kk = k - 1
+ do i=1,lonr
+ ii = kk*lonrbm + i
+ if (ii > lonr) ii = ii - lonr
+ rannum_tank(i,nrc+kk*maxrs,j) = rannum(ii+jj)
+ enddo
+ enddo
+ enddo
+ enddo
+ endif
+ endif
+!
+ if (me .eq. 0) then
+ print *,' seed0=',seed0,' idate=',idate,' wrk=',wrk
+ if (num_p3d .eq. 3) print *,' USING Ferrier-MICROPHYSICS'
+ if (num_p3d .eq. 4) print *,' USING ZHAO-MICROPHYSICS'
+ endif
+ if (fhour .eq. 0.0) then
+ do j=1,lats_node_r
+ do i=1,lonr
+! phy_f2d(i,j,num_p2d) = 0.0
+ phy_f2d(i,num_p2d,j) = 0.0
+ enddo
+ enddo
+ endif
+
+ if (ras) call ras_init(levs, me)
+
+ first = .false.
+!----------------------
+ endif ! if (first) done
+!----------------------
+!
+ dtphys = 3600.
+ nsphys = max(int((tstep+tstep)/dtphys+0.9999),1)
+ dtp = (tstep+tstep)/nsphys
+ dtf = 0.5*dtp
+!
+ if(lsfwd) dtf = dtp
+!
+ solhr = mod(phour+idate(1),cons_24)
+
+! ************** Ken Campana Stuff ********************************
+!... set switch for saving convective clouds
+ if(lscca.and.lsswr) then
+ clstp = 1100+min(fhswr,fhour,cons_99) !initialize,accumulate,convert
+ elseif(lscca) then
+ clstp = 0100+min(fhswr,fhour,cons_99) !accumulate,convert
+ elseif(lsswr) then
+ clstp = 1100 !initialize,accumulate
+ else
+ clstp = 0100 !accumulate
+ endif
+! ************** Ken Campana Stuff ********************************
+!
+!
+ iseed = mod(100.0*sqrt(fhour*3600),cons_1p0d9) + 1 + seed0
+
+
+ if (.not. newsas) then ! random number needed for RAS and old SAS
+ if (random_clds) then
+ call random_setseed(iseed)
+ call random_number(wrk2)
+ do nrc=1,nrcm
+ indxr(nrc) = max(1, min(nint(wrk2(nrc)*maxran)+1,maxran))
+ enddo
+ endif
+ endif
+!
+! doing ozone i/o and latitudinal interpolation to local gauss lats
+! ifozphys=.true.
+
+ if (ntoz .gt. 0) then
+ call ozinterpol(me,lats_node_r,lats_node_r,idate,fhour,
+ & jindx1,jindx2,ozplin,ozplout,ddy)
+
+!Moor call ozinterpol(lats_node_r,lats_node_r,idate,fhour,
+! & jindx1,jindx2,ozplin,ozplout,ddy,
+! & global_lats_r,lonsperlar)
+ endif
+
+! if (me == 0) write(0,*)' after ozinterpol'
+!!
+ pwatg = 0.
+ ptotg = 0.
+
+
+ do lan=1,lats_node_r
+ lat = global_lats_r(ipt_lats_node_r-1+lan)
+ lons_lat = lonsperlar(lat)
+ pwatp = 0.
+ rcs2_lan = rcs2_r(min(lat,latr-lat+1))
+ rcs_lan = sqrt(rcs2_lan)
+
+!$omp parallel do schedule(dynamic,1) private(lon)
+!$omp+private(sumq,xcp,hprime_v,swh_v,hlw_v,stc_v,smc_v,slc_v)
+!$omp+private(nlons_v,sinlat_v,coslat_v,ozplout_v,rannum_v)
+!$omp+private(prslk,prsl,prsik,prsi,phil,phii,dpshc,work1,tem)
+!$omp+private(gu,gv1,gd,gq,gphi,glam,gt,gtv,gr,vvel,gtvx,gtvy)
+!$omp+private(adt,adr,adu,adv,pgr,ugrd,vgrd,rqtk)
+!!$omp+private(adt,adr,adu,adv,pgr,rcs_v,ugrd,vgrd,rqtk)
+!$omp+private(phy_f3dv,phy_f2dv)
+!$omp+private(dt3dt_v,dq3dt_v,du3dt_v,dv3dt_v,upd_mf_v,dwn_mf_v)
+!$omp+private(det_mf_v,dkh_v,rnp_v)
+!$omp+private(njeff,item,jtem,ktem,i,j,k,n,kss)
+!!!$omp+private(temlon,temlat,lprnt,ipt)
+ do lon=1,lons_lat,ngptc
+!!
+ njeff = min(ngptc,lons_lat-lon+1)
+!!
+! lprnt = .false.
+!
+
+-------------------------------
+!---from MPAS
+ do k = 1, LEVS
+ do j = 1, njeff
+ jtem = lon-1+j
+ gu (j,k) =
+ gv1(j,k) =
+ gd (j,k) =
+ gt (j,k) = gtv(j,k) / (1.0 + fv*max(gr(j,k,1),qmin))
+ enddo
+ enddo
+!
+ do i=1,njeff
+ item = lon+i-1
+ gq(i) = for_gr_r_2(item,ksq,lan)
+ gphi(i) = for_gr_r_2(item,kspphi,lan)
+ glam(i) = for_gr_r_2(item,ksplam,lan)
+ enddo
+! Tracers
+ do n=1,ntrac
+ do k=1,levs
+ item = KSR-1+k+(n-1)*levs
+ do j=1,njeff
+ gr(j,k,n) = for_gr_r_2(lon-1+j,item,lan)
+ enddo
+ enddo
+ enddo
+-------------------------------
+
+!
+!
+ do i=1,njeff
+ phil(i,levs) = 0.0 ! forces calculation of geopotential in gbphys
+ pgr(i) = gq(i) * 1000.0 ! Convert from kPa to Pa for physics
+ prsi(i,1) = pgr(i)
+ dpshc(i) = 0.3 * prsi(i,1)
+!
+ nlons_v(i) = lons_lat
+ sinlat_v(i) = sinlat_r(lat)
+ coslat_v(i) = coslat_r(lat)
+! rcs_v(i) = sqrt(rcs2_lan)
+! rcs_v(i) = sqrt(rcs2_r(min(lat,latr-lat+1)))
+ enddo
+ do k=1,levs
+ do i=1,njeff
+ ugrd(i,k) = gu(i,k) * rcs_lan
+ vgrd(i,k) = gv1(i,k) * rcs_lan
+! ugrd(i,k) = gu(i,k) * rcs_v(i)
+! vgrd(i,k) = gv1(i,k) * rcs_v(i)
+ prsl(i,k) = prsl(i,k) * 1000.0
+ prsi(i,k+1) = prsi(i,k+1) * 1000.0
+ vvel(i,k) = vvel(i,k) * 1000.0 ! Convert from Cb/s to Pa/s
+ enddo
+ enddo
+
+!??????????????????
+ if (gen_coord_hybrid .and. thermodyn_id == 3) then
+ do i=1,ngptc
+ prslk(i,1) = 0.0 ! forces calculation of geopotential in gbphys
+ prsik(i,1) = 0.0 ! forces calculation of geopotential in gbphys
+ enddo
+ endif
+!
+ if (ntoz .gt. 0) then
+ do j=1,pl_coeff
+ do k=1,levozp
+ do i=1,ngptc
+!Moorthi ozplout_v(i,k,j) = ozplout(lon+i-1,k,j,lan)
+ ozplout_v(i,k,j) = ozplout(k,lan,j)
+ enddo
+ enddo
+ enddo
+ endif
+ do k=1,lsoil
+ do i=1,njeff
+ item = lon+i-1
+ smc_v(i,k) = sfc_fld%smc(item,k,lan)
+ stc_v(i,k) = sfc_fld%stc(item,k,lan)
+ slc_v(i,k) = sfc_fld%slc(item,k,lan)
+ enddo
+ enddo
+ do k=1,nmtvr
+ do i=1,njeff
+ hprime_v(i,k) = hprime(lon+i-1,k,lan)
+ enddo
+ enddo
+!!
+ do j=1,num_p3d
+ do k=1,levs
+ do i=1,njeff
+ phy_f3dv(i,k,j) = phy_f3d(lon+i-1,k,j,lan)
+ enddo
+ enddo
+ enddo
+ do j=1,num_p2d
+ do i=1,njeff
+ phy_f2dv(i,j) = phy_f2d(lon+i-1,j,lan)
+ enddo
+ enddo
+ if (.not. newsas) then
+ if (random_clds) then
+ do j=1,nrcm
+ do i=1,njeff
+ rannum_v(i,j) = rannum_tank(lon+i-1,indxr(j),lan)
+ enddo
+ enddo
+ else
+ do j=1,nrcm
+ do i=1,njeff
+ rannum_v(i,j) = 0.6 ! This is useful for debugging
+ enddo
+ enddo
+ endif
+ endif
+ do k=1,levs
+ do i=1,njeff
+ item = lon+i-1
+ swh_v(i,k) = swh(item,k,lan)
+ hlw_v(i,k) = hlw(item,k,lan)
+ enddo
+ enddo
+ if (ldiag3d) then
+ do n=1,6
+ do k=1,levs
+ do i=1,njeff
+ dt3dt_v(i,k,n) = dt3dt(lon+i-1,k,n,lan)
+ enddo
+ enddo
+ enddo
+ do n=1,4
+ do k=1,levs
+ do i=1,njeff
+ du3dt_v(i,k,n) = du3dt(lon+i-1,k,n,lan)
+ dv3dt_v(i,k,n) = dv3dt(lon+i-1,k,n,lan)
+ enddo
+ enddo
+ enddo
+ endif
+ if (ldiag3d .or. lggfs3d) then
+ do n=1,5+pl_coeff
+ do k=1,levs
+ do i=1,njeff
+ dq3dt_v(i,k,n) = dq3dt(lon+i-1,k,n,lan)
+ enddo
+ enddo
+ enddo
+ do k=1,levs
+ do i=1,njeff
+ upd_mf_v(i,k) = upd_mf(lon+i-1,k,lan)
+ dwn_mf_v(i,k) = dwn_mf(lon+i-1,k,lan)
+ det_mf_v(i,k) = det_mf(lon+i-1,k,lan)
+ enddo
+ enddo
+ endif
+ if (lggfs3d) then
+ do k=1,levs
+ do i=1,njeff
+ dkh_v(i,k) = dkh(lon+i-1,k,lan)
+ rnp_v(i,k) = rnp(lon+i-1,k,lan)
+ enddo
+ enddo
+ endif
+!
+! write(0,*)' calling gbphys kdt=',kdt,' lon=',lon,' lan=',lan
+! &,' nlons_v=',ntoz,ntcw,nmtvr,lonr,latr,jcap,ras
+! &,' tisfc=',sfc_fld%tisfc(lon,lan)
+! print *,' temp=',for_gr_r_2(lon,kst,lan)
+!
+ call gbphys &
+! --- inputs:
+ & ( njeff,ngptc,levs,lsoil,lsm,ntrac,ncld,ntoz,ntcw, &
+ & nmtvr,nrcm,levozp,lonr,latr,jcap,num_p3d,num_p2d, &
+ & kdt,lat,me,pl_coeff,nlons_v,ncw,flgmin,crtrh,cdmbgwd, &
+ & ccwf,dlqf,ctei_rm,clstp,dtp,dtf,fhour,solhr, &
+ & slag,sdec,cdec,sinlat_v,coslat_v,pgr,ugrd,vgrd, &
+ & gt,gr,vvel,prsi,prsl,prslk,prsik,phii,phil, &
+ & rannum_v,ozplout_v,pl_pres,dpshc, &
+ & hprime_v, xlon(lon,lan),xlat(lon,lan), &
+ & sfc_fld%slope (lon,lan), sfc_fld%shdmin(lon,lan), &
+ & sfc_fld%shdmax(lon,lan), sfc_fld%snoalb(lon,lan), &
+ & sfc_fld%tg3 (lon,lan), sfc_fld%slmsk (lon,lan), &
+ & sfc_fld%vfrac (lon,lan), sfc_fld%vtype (lon,lan), &
+ & sfc_fld%stype (lon,lan), sfc_fld%uustar(lon,lan), &
+ & sfc_fld%oro (lon,lan), flx_fld%coszen(lon,lan), &
+ & flx_fld%sfcdsw(lon,lan), flx_fld%sfcnsw(lon,lan), &
+ & flx_fld%sfcdlw(lon,lan), flx_fld%tsflw (lon,lan), &
+ & flx_fld%sfcemis(lon,lan), sfalb(lon,lan), &
+ & swh_v, hlw_v, &
+ & ras,pre_rad,ldiag3d,lggfs3d,lssav,lssav_cc, &
+ & bkgd_vdif_m,bkgd_vdif_h,bkgd_vdif_s,psautco,prautco, evpco, &
+ & wminco, &
+ & flipv,old_monin,cnvgwd,shal_cnv,sashal,newsas,cal_pre, &
+ & mom4ice,mstrat,trans_trac,nst_fcst,moist_adj,fscav, &
+ & thermodyn_id, sfcpress_id, gen_coord_hybrid, &
+! --- input/outputs:
+ & sfc_fld%hice (lon,lan), sfc_fld%fice (lon,lan), &
+ & sfc_fld%tisfc (lon,lan), sfc_fld%tsea (lon,lan), &
+ & sfc_fld%tprcp (lon,lan), sfc_fld%cv (lon,lan), &
+ & sfc_fld%cvb (lon,lan), sfc_fld%cvt (lon,lan), &
+ & sfc_fld%srflag(lon,lan), sfc_fld%snwdph(lon,lan), &
+ & sfc_fld%sheleg(lon,lan), sfc_fld%sncovr(lon,lan), &
+ & sfc_fld%zorl (lon,lan), sfc_fld%canopy(lon,lan), &
+ & sfc_fld%ffmm (lon,lan), sfc_fld%ffhh (lon,lan), &
+ & sfc_fld%f10m (lon,lan), flx_fld%srunoff(lon,lan), &
+ & flx_fld%evbsa (lon,lan), flx_fld%evcwa (lon,lan), &
+ & flx_fld%snohfa(lon,lan), flx_fld%transa(lon,lan), &
+ & flx_fld%sbsnoa(lon,lan), flx_fld%snowca(lon,lan), &
+ & flx_fld%soilm (lon,lan), flx_fld%tmpmin(lon,lan), &
+ & flx_fld%tmpmax(lon,lan), flx_fld%dusfc (lon,lan), &
+ & flx_fld%dvsfc (lon,lan), flx_fld%dtsfc (lon,lan), &
+ & flx_fld%dqsfc (lon,lan), flx_fld%geshem(lon,lan), &
+ & flx_fld%gflux (lon,lan), flx_fld%dlwsfc(lon,lan), &
+ & flx_fld%ulwsfc(lon,lan), flx_fld%suntim(lon,lan), &
+ & flx_fld%runoff(lon,lan), flx_fld%ep (lon,lan), &
+ & flx_fld%cldwrk(lon,lan), flx_fld%dugwd (lon,lan), &
+ & flx_fld%dvgwd (lon,lan), flx_fld%psmean(lon,lan), &
+ & flx_fld%bengsh(lon,lan), flx_fld%spfhmin(lon,lan), &
+ & flx_fld%spfhmax(lon,lan), &
+ & dt3dt_v, dq3dt_v, du3dt_v, dv3dt_v, &
+ & acv(lon,lan), acvb(lon,lan), acvt(lon,lan), &
+ & slc_v, smc_v, stc_v, &
+ & upd_mf_v, dwn_mf_v, det_mf_v, dkh_v, rnp_v, &
+ & phy_f3dv, phy_f2dv, &
+ & DLWSFC_cc(lon,lan), ULWSFC_cc(lon,lan), &
+ & DTSFC_cc(lon,lan), SWSFC_cc(lon,lan), &
+ & DUSFC_cc(lon,lan), DVSFC_cc(lon,lan), &
+ & DQSFC_cc(lon,lan), PRECR_cc(lon,lan), &
+
+ & nst_fld%xt(lon,lan), nst_fld%xs(lon,lan), &
+ & nst_fld%xu(lon,lan), nst_fld%xv(lon,lan), &
+ & nst_fld%xz(lon,lan), nst_fld%zm(lon,lan), &
+ & nst_fld%xtts(lon,lan), nst_fld%xzts(lon,lan), &
+ & nst_fld%d_conv(lon,lan), nst_fld%ifd(lon,lan), &
+ & nst_fld%dt_cool(lon,lan), nst_fld%Qrain(lon,lan), &
+! --- outputs:
+ & adt, adr, adu, adv, &
+ & sfc_fld%t2m (lon,lan), sfc_fld%q2m (lon,lan), &
+ & flx_fld%u10m (lon,lan), flx_fld%v10m (lon,lan), &
+ & flx_fld%zlvl (lon,lan), flx_fld%psurf (lon,lan), &
+ & flx_fld%hpbl (lon,lan), flx_fld%pwat (lon,lan), &
+ & flx_fld%t1 (lon,lan), flx_fld%q1 (lon,lan), &
+ & flx_fld%u1 (lon,lan), flx_fld%v1 (lon,lan), &
+ & flx_fld%chh (lon,lan), flx_fld%cmm (lon,lan), &
+ & flx_fld%dlwsfci(lon,lan), flx_fld%ulwsfci(lon,lan), &
+ & flx_fld%dswsfci(lon,lan), flx_fld%uswsfci(lon,lan), &
+ & flx_fld%dtsfci(lon,lan), flx_fld%dqsfci(lon,lan), &
+ & flx_fld%gfluxi(lon,lan), flx_fld%epi (lon,lan), &
+ & flx_fld%smcwlt2(lon,lan), flx_fld%smcref2(lon,lan), &
+!hchuang code change [+3L] 11/12/2007 : add 2D
+ & flx_fld%gsoil(lon,lan), flx_fld%gtmp2m(lon,lan), &
+ & flx_fld%gustar(lon,lan), flx_fld%gpblh(lon,lan), &
+ & flx_fld%gu10m(lon,lan), flx_fld%gv10m(lon,lan), &
+ & flx_fld%gzorl(lon,lan), flx_fld%goro(lon,lan), &
+
+ & XMU_cc(lon,lan), DLW_cc(lon,lan), DSW_cc(lon,lan), &
+ & SNW_cc(lon,lan), LPREC_cc(lon,lan), &
+
+ & nst_fld%Tref(lon,lan), nst_fld%z_c(lon,lan), &
+ & nst_fld%c_0 (lon,lan), nst_fld%c_d(lon,lan), &
+ & nst_fld%w_0 (lon,lan), nst_fld%w_d(lon,lan), &
+ & rqtk &! rqtkD
+! & bak_gr_r_2(lon,kap,lan), &! rqtkD
+ & )
+!!
+!!
+ prsi = prsi * 0.001 ! Convert from Pa to kPa
+ do k=1,lsoil
+ do i=1,njeff
+ item = lon + i - 1
+ sfc_fld%smc(item,k,lan) = smc_v(i,k)
+ sfc_fld%stc(item,k,lan) = stc_v(i,k)
+ sfc_fld%slc(item,k,lan) = slc_v(i,k)
+ enddo
+ enddo
+!!
+ do j=1,num_p3d
+ do k=1,levs
+ do i=1,njeff
+ phy_f3d(lon+i-1,k,j,lan) = phy_f3dv(i,k,j)
+ enddo
+ enddo
+ enddo
+ do j=1,num_p2d
+ do i=1,njeff
+ phy_f2d(lon+i-1,j,lan) = phy_f2dv(i,j)
+ enddo
+ enddo
+!
+ if (ldiag3d) then
+ do n=1,6
+ do k=1,levs
+ do i=1,njeff
+ dt3dt(lon+i-1,k,n,lan) = dt3dt_v(i,k,n)
+ enddo
+ enddo
+ enddo
+ do n=1,4
+ do k=1,levs
+ do i=1,njeff
+ du3dt(lon+i-1,k,n,lan) = du3dt_v(i,k,n)
+ dv3dt(lon+i-1,k,n,lan) = dv3dt_v(i,k,n)
+ enddo
+ enddo
+ enddo
+ endif
+ if (ldiag3d .or. lggfs3d) then
+ do n=1,5+pl_coeff
+ do k=1,levs
+ do i=1,njeff
+ dq3dt(lon+i-1,k,n,lan) = dq3dt_v(i,k,n)
+ enddo
+ enddo
+ enddo
+ do k=1,levs
+ do i=1,njeff
+ upd_mf(lon+i-1,k,lan) = upd_mf_v(i,k)
+ dwn_mf(lon+i-1,k,lan) = dwn_mf_v(i,k)
+ det_mf(lon+i-1,k,lan) = det_mf_v(i,k)
+ enddo
+ enddo
+ endif
+ if (lggfs3d) then
+ do k=1,levs
+ do i=1,njeff
+ dkh(lon+i-1,k,lan) = dkh_v(i,k)
+ rnp(lon+i-1,k,lan) = rnp_v(i,k)
+ enddo
+ enddo
+ endif
+!
+!---------------------------
+ enddo ! lon loop
+!---------------------------
+!
+!
+ ptotj(lan) = 0.
+ do j=1,lons_lat
+ ptotj(lan) = ptotj(lan) + gq_save(j,lan)
+ pwatp = pwatp + flx_fld%pwat(j,lan)
+! print *,' kdt=',kdt,' pwatp=',pwatp,' pwat=',flx_fld%pwat(j,lan)
+! &,' j=',j
+ enddo
+ pwatj(lan) = pwatp*grav/(2.*lonsperlar(lat)*1.e3)
+ ptotj(lan) = ptotj(lan)/(2.*lonsperlar(lat))
+
+!---------------------------
+ enddo ! lan loop
+!---------------------------
+!!
+!
+!!?????????????????????????????????????????????????
+ call excha(lats_nodes_r,global_lats_r,ptotj,pwatj,ptotg,pwatg)
+ sumwa = 0.
+ sumto = 0.
+ do lat=1,latr
+ sumto = sumto + wgt_r(min(lat,latr-lat+1))*ptotg(lat)
+ sumwa = sumwa + wgt_r(min(lat,latr-lat+1))*pwatg(lat)
+! print *,' kdt=',kdt,' lat=',lat,' sumwa=',sumwa,' sumto=',sumto,
+! &' ptotg=',ptotg(lat),' pwatg=',pwatg(lat)
+ enddo
+cjfe
+cjfe write(70+me,*) sumto,sumwa,kdt
+ pdryg = sumto - sumwa
+!!
+ if(pdryini == 0.) pdryini = pdryg
+
+ if( gen_coord_hybrid ) then ! hmhj
+ pcorr = (pdryini-pdryg) * sqrt(2.) ! hmhj
+ else ! hmhj
+ pcorr = (pdryini-pdryg) / sumto * sqrt(2.)
+ endif ! hmhj
+!!?????????????????????????????????????????????????
+
+ return
+ end
Modified: branches/atmos_physics_gfs/src/core_atmos_physics_gfs/physics_gfs/gloopr.f
===================================================================
--- branches/atmos_physics_gfs/src/core_atmos_physics_gfs/physics_gfs/gloopr.f 2012-05-09 22:12:47 UTC (rev 1884)
+++ branches/atmos_physics_gfs/src/core_atmos_physics_gfs/physics_gfs/gloopr.f 2012-05-10 01:07:10 UTC (rev 1885)
@@ -4,7 +4,7 @@
subroutine gloopr
!---input
- & (lonsperlar,global_lats_r,phour,xlon,xlat,kdt,
+ & (phour,kdt,lonsperlar,global_lats_r,xlon,xlat,
& slmsk,sheleg, zorl, tsea,
& alvsf, alnsf, alvwf, alnwf, facsf, facwf,
& cv, cvt, cvb, fice, tisfc, sncovr, snoalb,
@@ -57,13 +57,13 @@
& hprime(lonr,nmtvr,lats_node_r), phour, &
& phy_f3d(lonr,levs,num_p3d,lats_node_r)
-! --prsi : model level pressure in Pasca
-! --prsl : model layer mean pressure in Pasca
+! --prsi : model level pressure in centipar
+! --prsl : model layer mean pressure in centibar
! --gt : model layer mean temperature in k
! --gr : layer specific humidity in gm/gm
! --gr1 : layer tracer (ozne and cloud water) mass mixing ratio
-! --vvel : layer mean vertical velocity in Pascal/sec
- real (kind=kind_phys), intent(inout) :: &
+! --vvel : layer mean vertical velocity in centibar/sec
+ real (kind=kind_phys), intent(in) ::
& prsi(lonr,levp1,lats_node_r),prsl(lonr,levs,lats_node_r), &
& gt(lonr,levs,lats_node_r),gr(lonr,levs,lats_node_r), &
& gr1(lonr,levs,ntrac-1,lats_node_r),vvel(lonr,levs,lats_node_r)
@@ -354,18 +354,19 @@
lat = global_lats_r(ipt_lats_node_r-1+j)
lons_lat = lonsperlar(lan)
-! --- vertical structure variable prslk and minimum water vapor mixing ratio
+! --- vertical structure variable prslk
+! and minimum water vapor mixing ratio
do k=1,levr
do i=1,lons_lat
- prslk(i,k,lan) = prslk(i,k,lan)**con_rocp
- gr(i,k,lan) = max(qmin,gr(i,k,lan))
+ prslk(i,k,lan) = (0.01*prsl(i,k,lan))**con_rocp
+!! gr(i,k,lan) = max(qmin,gr(i,k,lan))
enddo
enddo
!!
!$omp parallel do schedule(dynamic,1) private(lon,i,j,k)
!$omp+private(cldcov_v,fluxr_v,f_ice,f_rain,r_rime)
-!$omp+private(prslk,flgmin_v,hlw_v,swh_v)
+!$omp+private(flgmin_v,hlw_v,swh_v)
!$omp+private(njeff,n,item,jtem,ks,work1,work2)
!$omp+private(icsdsw,icsdlw)
!$omp+private(lprnt,ipt)
@@ -374,6 +375,7 @@
!--------------------
NJEFF = MIN(ngptc,lons_lat-lon+1)
lprnt = .false.
+ ipt=lon !diagnostic printout point
!
do k=1,nfxr
do j=1,njeff
</font>
</pre>