!at tune step========================================================= ! ! description: ! ! ! ! gbphys is the driver subroutine to invoke GFS AM physics ! ! (except radiation but radiative heating is applied here) ! ! at physics time steps ! ! ! ! usage: ! ! ! ! call gbphys ! ! inputs: ! ! ( im,ix,levs,lsoil,lsm,ntrac,ncld,ntoz,ntcw,ntke, ! ! nmtvr,nrcm,ko3,lonr,latr,jcap, ! ! num_p3d,num_p2d,npdf3d,ncnvcld3d, ! ! kdt,lat,me,pl_coeff,nlons,ncw,flgmin,crtrh,cdmbgwd, ! ! ccwf,dlqf,ctei_rm,clstp,cgwf,prslrd0,ral_ts,dtp,dtf,fhour, ! ! solhr,slag,sdec,cdec,sinlat,coslat,pgr,ugrs,vgrs, ! ! tgrs,qgrs,vvel,prsi,prsl,prslk,prsik,phii,phil, ! ! rann,prdout,poz,dpshc,hprime,xlon,xlat, ! ! slope,shdmin,shdmax,snoalb,tg3,slmsk,vfrac, ! ! vtype,stype,uustar,oro,oro_uf,coszen,sfcdsw,sfcnsw, ! ! sfcnirbmd,sfcnirdfd,sfcvisbmd,sfcvisdfd, ! ! sfcnirbmu,sfcnirdfu,sfcvisbmu,sfcvisdfu, ! ! slimskin_cpl,ulwsfcin_cpl, ! ! dusfcin_cpl,dvsfcin_cpl,dtsfcin_cpl,dqsfcin_cpl, ! ! sfcdlw,tsflw,sfcemis,sfalb,swh,swhc,hlw,hlwc,hlwd,lsidea, ! ! ras,pre_rad,ldiag3d,lgocart,lssav,lssav_cpl ! ! xkzm_m,xkzm_h,xkzm_s,psautco,prautco,evpco,wminco, ! ! pdfcld,shcnvcw,sup,redrag,hybedmf,dspheat, ! ! flipv,old_monin,cnvgwd,shal_cnv, ! ! imfshalcnv,imfdeepcnv,cal_pre, ! ! mom4ice,mstrat,trans_trac,nstf_name,moist_adj, ! ! thermodyn_id, sfcpress_id, gen_coord_hybrid,levr,adjtrc,nnp,! ! cscnv,nctp,do_shoc,shocaftcnv,ntot3d,ntot2d, ! ! input/outputs: ! ! hice,fice,tisfc,tsea,tprcp,cv,cvb,cvt, ! ! srflag,snwdph,weasd,sncovr,zorl,canopy, ! ! ffmm,ffhh,f10m,srunoff,evbsa,evcwa,snohfa, ! ! transa,sbsnoa,snowca,soilm,tmpmin,tmpmax, ! ! dusfc,dvsfc,dtsfc,dqsfc,totprcp,gflux, ! ! dlwsfc,ulwsfc,suntim,runoff,ep,cldwrk, ! ! dugwd,dvgwd,psmean,cnvprcp,spfhmin,spfhmax,rain,rainc, ! ! dt3dt,dq3dt,du3dt,dv3dt,dqdt_v,cnvqc_v,acv,acvb,acvt, ! ! slc,smc,stc,upd_mf,dwn_mf,det_mf,phy_f3d,phy_f2d, ! ! dusfc_cpl, dvsfc_cpl, dtsfc_cpl, dqsfc_cpl, ! ! dlwsfc_cpl,dswsfc_cpl,dnirbm_cpl,dnirdf_cpl, ! ! dvisbm_cpl,dvisdf_cpl,rain_cpl, nlwsfc_cpl,nswsfc_cpl, ! ! nnirbm_cpl,nnirdf_cpl,nvisbm_cpl,nvisdf_cpl,snow_cpl, ! ! xt,xs,xu,xv,xz,zm,xtts,xzts,d_conv,ifd,dt_cool,Qrain, ! ! phy_fctd, ! ! outputs: ! ! gt0,gq0,gu0,gv0,t2m,q2m,u10m,v10m, ! ! zlvl,psurf,hpbl,pwat,t1,q1,u1,v1, ! ! chh,cmm,dlwsfci,ulwsfci,dswsfci,uswsfci,dusfci,dvsfci, ! ! dtsfci,dqsfci,gfluxi,epi,smcwlt2,smcref2,wet1, ! ! dusfci_cpl,dvsfci_cpl,dtsfci_cpl,dqsfci_cpl, ! ! dlwsfci_cpl,dswsfci_cpl, ! ! dnirbmi_cpl,dnirdfi_cpl,dvisbmi_cpl,dvisdfi_cpl, ! ! nlwsfci_cpl,nswsfci_cpl, ! ! nnirbmi_cpl,nnirdfi_cpl,nvisbmi_cpl,nvisdfi_cpl, ! ! t2mi_cpl,q2mi_cpl, ! ! u10mi_cpl,v10mi_cpl,tseai_cpl,psurfi_cpl, ! ! tref, z_c, c_0, c_d, w_0, w_d, rqtk, ! ! hlwd,lsidea ) ! ! ! ! subprograms called: ! ! ! ! get_prs, dcyc2t2_pre_rad (testing), dcyc2t3, sfc_diff, ! ! sfc_ocean,sfc_drv, sfc_land, sfc_sice, sfc_diag, moninp1, ! ! moninp, moninq1, moninq, gwdps, ozphys, get_phi, ! ! sascnv, sascnvn, rascnv, cs_convr, gwdc, shalcvt3,shalcv,! ! shalcnv, cnvc90, lrgscl, gsmdrive, gscond, precpd, ! ! progt2. ! ! ! ! ! ! program history log: ! ! 19xx - ncep mrf/gfs ! ! 2002 - s. moorthi modify and restructure and add Ferrier ! ! microphysics as an option ! ! 200x - h. juang modify (what?) ! ! nov 2004 - x. wu modify sea-ice model ! ! may 2005 - s. moorthi modify and restructure ! ! 2005 - s. lu modify to include noah lsm ! ! oct 2006 - h. wei modify lsm options to include both ! ! noah and osu lsms. ! ! 2006 - s. moorthi added a. johansson's convective gravity ! ! wave parameterization code ! ! 2007 - s. moorthi added j. han's modified pbl/sas options ! ! dec 2007 - xu li modified the operational version for ! ! nst model ! ! 2008 - s. moorthi applied xu li's nst model to new gfs ! ! mar 2008 - y.-t. hou added sunshine duration var (suntim) as ! ! an input/output argument. ! ! 2008 - jun wang added spfhmax/spfhmin as input/output. ! ! apr 2008 - y.-t. hou added lw sfc emissivity var (sfcemis), ! ! define the lw sfc dn/up fluxes in two forms: atmos! ! and ground. also changed sw sfc net flux direction! ! (positive) from ground -> atmos to the direction ! ! of atmos -> ground. recode the program and add ! ! program documentation block. ! 2008/ - s. moorthi and y.t. hou upgraded the code to more ! ! 2009 modern form and changed all the inputs to MKS units.! ! feb 2009 - s. moorthi upgraded to add Hochun's gocart changes ! ! jul 2009 - s. moorthi added rqtk for sela's semi-lagrangian ! ! aug 2009 - s. moorthi added j. han and h. pan updated shallow ! ! convection package ! ! sep 2009 - s. moorthi updated for the mcica (rrtm3) radiation ! ! dec 2010 - sarah lu lgocart added to input arg; ! ! compute dqdt_v if inline gocart is on ! ! feb 2011 - sarah lu add the option to update surface diag ! ! fields (t2m,q2m,u10m,v10m) at the end ! ! Jun 2011 - s. moorthi and Xu Li - updated the nst model ! ! ! ! sep 2011 - sarah lu correct dqdt_v calculations ! ! apr 2012 - henry juang add idea ! ! sep 2012 - s. moorthi merge with operational version ! ! Mar 2013 - Jun Wang set idea heating rate to tmp tendency ! ! May 2013 - Jun Wang tmp updated after idea phys ! ! Jun 2013 - s. moorthi corrected a bug in 3d diagnostics for T ! ! Aug 2013 - s. moorthi updating J. Whitekar's changes related ! ! to stochastic physics perturnbation ! ! Oct 2013 - Xingren Wu add dusfci/dvsfci ! ! Mar 2014 - Xingren Wu add "_cpl" for coupling ! ! Mar 2014 - Xingren Wu add "nir/vis beam and nir/vis diff" ! ! Apr 2014 - Xingren Wu add "NET LW/SW including nir/vis" ! ! Jan 2014 - Jun Wang merge Moorthi's gwdc change and H.Juang ! ! and F. Yang's energy conversion from GWD! ! jan 2014 - y-t hou revised sw sfc spectral component fluxes! ! for coupled mdl, added estimation of ocean albedo ! ! without ice contamination. ! ! Jun 2014 - Xingren Wu update net SW fluxes over the ocean ! ! (no ice contamination) ! ! Jul 2014 - Xingren Wu add Sea/Land/Ice Mask - slmsk_cpl ! ! Jul 2014 - s. moorthi merge with standalone gfs and cleanup ! ! Aug 2014 - s. moorthi add tracer fixer ! ! Sep 2014 - Sarah Lu disable the option to compute tracer ! ! scavenging in GFS phys (set fscav=0.) ! ! Dec 2014 - Jun Wang add cnvqc_v for gocart ! ! ==================== defination of variables ==================== ! ! --- 2014 - D. Dazlich Added Chikira-Sugiyama (CS) convection ! ! as an option in opr GFS. ! ! Apr 2015 S. Moorthi Added CS scheme to NEMS/GSM ! ! Jun 2015 S. Moorthi Added SHOC to NEMS/GSM ! ! Aug 2015 - Xu Li change nst_fcst to be nstf_name ! ! and introduce depth mean SST ! ! Sep 2015 - Xingren Wu remove oro_cpl & slmsk_cpl ! ! Sep 2015 - Xingren Wu add sfc_cice ! ! Sep 2015 - Xingren Wu connect CICE output to sfc_cice ! ! Jan 2016 - P. Tripp NUOPC/GSM merge ! ! Mar 2016 - J. Han - add ncnvcld3d integer ! for convective cloudiness enhancement ! Mar 2016 - J. Han - change newsas & sashal to imfdeepcnv ! & imfshalcnv, respectively ! Mar 2016 F. Yang add pgr to rayleigh damping call ! ! Mar 2016 S. Moorthi add ral_ts ! ! ==================== definition of variables ==================== ! ! ! ! inputs: size ! ! ix, im - integer, horiz dimention and num of used pts 1 ! ! levs - integer, vertical layer dimension 1 ! ! lsoil - integer, number of soil layers 1 ! ! lsm - integer, flag for land surface model to use 1 ! ! =0 for osu lsm; =1 for noah lsm ! ! ntrac - integer, number of tracers 1 ! ! ncld - integer, number of cloud species 1 ! ! ntoz - integer, ozone location in the tracer array 1 ! ! ntcw - integer, cloud condensate location in the tracer 1 ! ! array 1 ! ! ntke - integer, tke location in the tracer array 1 ! ! nmtvr - integer, number of topographic variables such as 1 ! ! variance etc used in the GWD parameterization ! ! nrcm - integer, second dimension for the random number 1 ! ! array rann ! ! ko3 - integer, number of layers for ozone data 1 ! ! lonr,latr- integer, number of lon/lat points 1 ! ! jcap - integer, number of spectral wave trancation 1 ! ! used only by sascnv shalcnv ! ! num_p3d - integer, number of 3D arrays needed for 1 ! ! microphysics ! ! num_p2d - integer, number of 2D arrays needed for 1 ! ! microphysics ! ! npdf3d - integer, number of 3d arrays associated with pdf 1 ! ! based clouds/microphysics ! ! ncnvcld3d- integer, number of 3d arrays associated with 1 ! ! convective cloudiness enhancement ! ! kdt -integer, number of the current time step 1 ! ! lat -integer, latitude index - used for debug prints 1 ! ! me -integer, pe number - used for debug prints 1 ! ! pl_coeff - integer, number coefficients in ozone forcing 1 ! ! nlons - integer, number of total grid points in a latitude ! ! circle through a point im ! ! ncw - integer, range of droplet number concentrations for ! ! Ferrier microphysics 2 ! ! flgmin - real, range of minimum large ice fraction for ! ! Ferrier microphys 2 ! ! crtrh - real, critical relative humidity at the surface, PBL ! ! top and at the top of the atmosphere 3 ! ! cdmbgwd - real, multiplication factors for cdmb and gwd 2 ! ! ccwf - real, multiplication factor for critical cloud ! ! workfunction for RAS 2 ! ! dlqf - real, factor for cloud condensate detrainment from ! ! cloud edges (RAS) 2 ! ! ctei_rm - real, critical cloud top entrainment instability 2 ! ! criteria (used if mstrat=.true.) ! ! clstp - real, index used by cnvc90 (for convective clouds)1 ! ! legacy stuff - does not affect forecast ! ! cgwf - real, multiplication factor for convective GWD 2 ! ! prslrd0 - real, pressure level (Pa) from which Rayleigh Damping ! ! is applied 1 ! ! ral_ts - real time scale for Rayleigh damping in days 1 ! ! dtp - real, physics time step in seconds 1 ! ! dtf - real, dynamics time step in seconds 1 ! ! fhour - real, forecast hour 1 ! ! solhr - real, fcst hour at the end of prev time step 1 ! ! slag - real, equation of time ( radian ) 1 ! ! sdec,cdec- real, sin and cos of the solar declination angle 1 ! ! sinlat - real, sin of latitude im ! ! coslat - real, cos of latitude im ! ! pgr - real, surface pressure (Pa) im ! ! ugrs,vgrs- real, u/v component of layer wind ix,levs ! ! tgrs - real, layer mean temperature ( k ) ix,levs ! ! qgrs - real, layer mean tracer concentration ix,levs,ntrac! ! vvel - real, layer mean vertical velocity (Pa/s) ix,levs ! ! prsi - real, pressure at layer interfaces ix,levs+1 ! prsl - real, mean layer presure ix,levs ! ! prsik - real, Exner function at layer interface ix,levs+1 ! prslk - real, Exner function at layer ix,levs ! ! phii - real, interface geopotential (m^2/s^2) ix,levs+1! ! phil - real, layer geopotential (m^2/s^2) ix,levs ! ! rann - real, random number array (0-1) ix,nrcm ! ! prdout - real, ozone forcing data ix,ko3,pl_coeff! ! poz - real, ozone forcing data level pressure (ln(Pa)) ko3 ! ! dpshc - real, maximum pressure depth for shallow convection im ! ! hprime - real, orographic std dev ix,nmtvr! ! xlon,xlat- real, longitude and latitude ( radian ) im ! ! slope - real, sfc slope type for lsm im ! ! shdmin - real, min fractional coverage of green veg im ! ! shdmax - real, max fractnl cover of green veg (not used) im ! ! snoalb - real, max snow albedo over land (for deep snow) im ! ! tg3 - real, deep soil temperature im ! ! slmsk - real, sea/land/ice mask (=0/1/2) im ! ! vfrac - real, vegetation fraction im ! ! vtype - real, vegetation type im ! ! stype - real, soil type im ! ! uustar - real, boundary layer parameter im ! ! oro - real, orography im ! ! oro_uf - real, unfiltered orography im ! ! coszen - real, avg cosz over daytime sw radiation interval im ! ! sfcdsw - real, total sky sfc downward sw flux ( w/m**2 ) im ! ! sfcnsw - real, total sky sfc netsw flx into ground(w/m**2) im ! ! sfcdlw - real, total sky sfc downward lw flux ( w/m**2 ) im ! ! tsflw - real, sfc air (layer 1) temp over lw interval (k) im ! ! sfcemis - real, sfc lw emissivity ( fraction ) im ! ! sfalb - real, mean sfc diffused sw albedo im ! ! sfcnirbmu- real, sfc nir-beam sw upward flux (w/m2) im ! ! sfcnirdfu- real, sfc nir-diff sw upward flux (w/m2) im ! ! sfcvisbmu- real, sfc uv+vis-beam sw upward flux (w/m2) im ! ! sfcvisdfu- real, sfc uv+vis-diff sw upward flux (w/m2) im ! ! sfcnirbmd- real, sfc nir-beam sw downward flux (w/m2) im ! ! sfcnirdfd- real, sfc nir-diff sw downward flux (w/m2) im ! ! sfcvisbmd- real, sfc uv+vis-beam sw downward flux (w/m2) im ! ! sfcvisdfd- real, sfc uv+vis-diff sw downward flux (w/m2) im ! ! slimskin_cpl - real, im ! ! ulwsfcin_cpl - real, im ! ! dusfcin_cpl - real, im ! ! dvsfcin_cpl - real, im ! ! dtsfcin_cpl - real, im ! ! dqsfcin_cpl - real, im ! ! swh - real, total sky sw heating rates ( k/s ) ix,levs ! ! swhc - real, clear sky sw heating rates ( k/s ) ix,levs ! ! hlw - real, total sky lw heating rates ( k/s ) ix,levs ! ! hlwc - real, clear sky lw heating rates ( k/s ) ix,levs ! ! hlwd - real, idea sky lw heating rates ( k/s ) ix,levs ! ! ras - logical, flag for ras convection scheme 1 ! ! cscnv - logical, flag for Chikira-Sugiyama convection 1 ! ! nctp - integer, number of cloud types in CS scheme 1 ! ! do_shoc - logical, flag for SHOC 1 ! ! shocaftcnv - logical, flag for SHOC 1 ! ! ntot3d - integer, number of total 3d fields for phy_f3d 1 ! ! ntot2d - integer, number of total 2d fields for phy_f2d 1 ! ! pre_rad - logical, flag for testing purpose 1 ! ! ldiag3d - logical, flag for 3d diagnostic fields 1 ! ! lgocart - logical, flag for 3d diagnostic fields for gocart 1 ! ! lssav - logical, flag controls data store and output 1 ! ! lssav_cpl- logical, flag for save data for A/O/I coupling 1 ! ! flipv - logical, flag for vertical direction flip (ras) 1 ! ! xkzm_m - real, background vertical diffusion for momentum 1 ! ! xkzm_h - real, background vertical diffusion for heat, q 1 ! ! xkzm_s - real, sigma threshold for background mom. diffusn 1 ! ! psautco - real, auto conversion coeff from ice to snow 2 ! ! prautco - real, auto conversion coeff from cloud to rain 2 ! ! evpco - real, coeff for evaporation of largescale rain 1 ! ! wminco - real, water and ice minimum threshold for Zhao 1 ! ! pdfcld - logical, flag for pdfcld 1 ! ! shcnvcw - logical, flag for shallow convective cloud 1 ! ! sup - real, supsaturation for ho. nucleation of ice 1 ! ! redrag - logical, flag for reduced drag coeff. over sea 1 ! ! hybedmf - logical, flag for hybrid edmf pbl scheme 1 ! ! dspheat - logical, flag for tke dissipative heating 1 ! ! old_monin- logical, flag for diff monin schemes 1 ! ! cnvgwd - logical, flag for conv gravity wave drag 1 ! ! shal_cnv - logical, flag for calling shallow convection 1 ! ! imfshalcnv - integer, flag for mass-flux shallow conv scheme 1 ! ! 1: July 2010 version of mass-flux shallow conv scheme ! current operational version as of 2016 ! 2: scale- & aerosol-aware mass-flux shallow conv scheme (2017) ! 0: modified Tiedtke's eddy-diffusion shallow conv scheme ! -1: no shallow convection used ! imfdeepcnv - integer, flag for mass-flux deep conv scheme 1 ! ! 1: July 2010 version of SAS conv scheme ! current operational version as of 2016 ! 2: scale- & aerosol-aware mass-flux deep conv scheme (2017) ! 0: old SAS Convection scheme before July 2010 ! cal_pre - logical, flag controls precip type algorithm 1 ! ! mom4ice - logical, flag controls mom4 sea-ice 1 ! ! mstrat - logical, flag for moorthi approach for stratus 1 ! ! trans_trac-logical, flag for convective transport of tracers 1 ! ! nstf_name -integer array, NSST related flag parameters 1 ! ! nstf_name(1) : 0 = NSSTM off 1 ! ! 1 = NSSTM on but uncoupled 1 ! ! 2 = NSSTM on and coupled 1 ! ! nstf_name(2) : 1 = NSSTM spin up on 1 ! ! 0 = NSSTM spin up off 1 ! ! nstf_name(3) : 1 = NSST analysis on 1 ! ! 0 = NSSTM analysis off 1 ! ! nstf_name(4) : zsea1 in mm 1 ! ! nstf_name(5) : zsea2 in mm 1 ! ! moist_adj- logical, flag for moist convective adjustment 1 ! ! thermodyn_id - integer, valid for GFS only for get_prs/phi 1 ! ! sfcpress_id - integer, valid for GFS only for get_prs/phi 1 ! ! gen_coord_hybrid - logical for pk=ak+bk*ps+ck*theta (Henry) 1 ! ! levr - integer, the number of layers GFS Radiative heating calculated at 1 ! ! adjtrc - real, dynamics adjustments to tracers ntrac ! ! nnp - integer, physics substep number 1 ! ! ! ! input/outputs: ! ! hice - real, sea-ice thickness im ! ! fice - real, sea-ice concentration im ! ! tisfc - real, sea-ice temperature im ! ! tsea - real, ground surface temperature ( k ) im ! ! tprcp - real, total precipitation im ! ! the following three variables do not affect the forecast ! ! cv, -real, convective clouds amountt im ! ! cvb -real, convective clouds base pressure (kPa) im ! ! cvt -real, convective clouds top pressure (kPa) im ! ! srflag - real, snow/rain flag for precipitation im ! ! snwdph - real, actual snow depth (mm) over land/sea ice im ! ! weasd - real, water equiv of accumulated snow depth (kg/m**2) ! ! over land and sea ice im ! ! sncovr - real, snow cover over land im ! ! zorl - real, surface roughness im ! ! canopy - real, canopy water im ! ! ffmm - real, fm parameter from PBL scheme im ! ! ffhh - real, fh parameter from PBL scheme im ! ! f10m - real, fm at 10m im ! ! srunoff - real, surface water runoff (from lsm) im ! ! evbsa - real, noah lsm diagnostics im ! ! evcwa - real, noah lsm diagnostics im ! ! snohfa - real, noah lsm diagnostics im ! ! transa - real, noah lsm diagnostics im ! ! sbsnoa - real, noah lsm diagnostics im ! ! snowca - real, noah lsm diagnostics im ! ! soilm - real, soil moisture im ! ! tmpmin - real, min temperature at 2m height (k) im ! ! tmpmax - real, max temperature at 2m height (k) im ! ! dusfc - real, u component of surface stress im ! ! dvsfc - real, v component of surface stress im ! ! dtsfc - real, sensible heat flux (w/m2) im ! ! dqsfc - real, latent heat flux (w/m2) im ! ! totprcp - real, accumulated total precipitation (kg/m2) im ! ! gflux - real, groud conductive heat flux im ! ! dlwsfc - real, time accumulated sfc dn lw flux ( w/m**2 ) im ! ! ulwsfc - real, time accumulated sfc up lw flux ( w/m**2 ) im ! ! suntim - real, sunshine duration time (s) im ! ! runoff - real, total water runoff im ! ! ep - real, potential evaporation im ! ! cldwrk - real, cloud workfunction (valid only with sas) im ! ! dugwd - real, vertically integrated u change by OGWD im ! ! dvgwd - real, vertically integrated v change by OGWD im ! ! psmean - real, surface pressure (kPa) im ! ! cnvprcp - real, accumulated convective precipitation (kg/m2)im ! ! spfhmin - real, minimum specific humidity im ! ! spfhmax - real, maximum specific humidity im ! ! rain - real, total rain at this time step im ! ! rainc - real, convective rain at this time step im ! ! dt3dt - real, temperature change due to physics ix,levs,6 ! ! dq3dt - real, moisture change due to physics ix,levs,5+pl_coeff! ! du3dt - real, u momentum change due to physics ix,levs,4 ! ! dv3dt - real, v momentum change due to physics ix,levs,4 ! ! dqdt_v - real, total moisture tendency (kg/kg/s) ix,levs ! ! cnvqc_v - real, total convective conensate (kg/kg) ix,levs ! ! acv - real, array containing accumulated convective clouds im ! ! acvb,acvt - real, arrays used by cnvc90 im ! ! slc - real, liquid soil moisture ix,lsoil! ! smc - real, total soil moisture ix,lsoil! ! stc - real, soil temperature ix,lsoil! ! upd_mf - real, convective updraft mass flux ix,levs ! ! dwn_mf - real, convective downdraft mass flux ix,levs ! ! det_mf - real, convective detrainment mass flux ix,levs ! ! ------- not used below ----------- ! dkh - real, vertical diffusion coefficient (gocart) ix,levs ! ! rnp - real, n cloud precipitation rate (gocart) ix,levs ! ! ------- not used above ----------- ! phy_f3d - real, 3d arrays saved for restart ix,levs,num_p3d! ! +npdf3d! ! phy_f2d - real, 2d arrays save for restart ix,num_p2d! ! dusfc_cpl - real, sfc u-momentum flux for A/O/I coupling im ! ! dvsfc_cpl - real, sfc v-momentum flux for A/O/I coupling im ! ! dtsfc_cpl - real, sfc sensible heat flux for A/O/I coupling im ! ! dqsfc_cpl - real, sfc latent heat flux for A/O/I coupling im ! ! dlwsfc_cpl- real, sfc dnwd lw flux (w/m**2) for A/O/I coupling im ! ! dswsfc_cpl- real, sfc dnwd sw flux (w/m**2) for A/O/I coupling im ! ! dnirbm_cpl- real, sfc nir beam dnwd sw rad flux (w/m**2) im ! ! dnirdf_cpl- real, sfc nir diff dnwd sw rad flux (w/m**2) im ! ! dvisbm_cpl- real, sfc uv+vis beam dnwd sw rad flux (w/m**2) im ! ! dvisdf_cpl- real, sfc uv+vis diff dnwd sw rad flux (w/m**2) im ! ! nlwsfc_cpl- real, net dnwd lw flux (w/m**2) for A/O/I coupling im ! ! nswsfc_cpl- real, net dnwd sw flux (w/m**2) for A/O/I coupling im ! ! nnirbm_cpl- real, net nir beam dnwd sw rad flux (w/m**2) im ! ! nnirdf_cpl- real, net nir diff dnwd sw rad flux (w/m**2) im ! ! nvisbm_cpl- real, net uv+vis beam dnwd sw rad flux (w/m**2) im ! ! nvisdf_cpl- real, net uv+vis diff dnwd sw rad flux (w/m**2) im ! ! rain_cpl - real, total precipitation rain for A/O/I coupling im ! ! snow_cpl - real, total precipitation snow for A/O/I coupling im ! ! ! xt - real, heat content in DTL im ! ! xs - real, salinity content in DTL im ! ! xu - real, u-current content in DTL im ! ! xv - real, v-current content in DTL im ! ! xz - real, DTL thickness im ! ! zm - real, MXL thickness im ! ! xtts - real, d(xt)/d(ts) im ! ! xzts - real, d(xz)/d(ts) im ! ! d_conv - real, thickness of Free Convection Layer (FCL) im ! ! ifd - real, index to start DTM run or not im ! ! dt_cool - real, Sub-layer cooling amount im ! ! Qrain - real, sensible heat flux due to rainfall (watts) im ! ! phy_fctd - real, cloud base mass flux for CScnv ix,nctp ! ! ! ! outputs: ! ! gt0 - real, updated temperature ix,levs ! ! gq0 - real, updated tracers ix,levs,ntrac! ! gu0 - real, updated zonal wind ix,levs ! ! gv0 - real, update meridional wind ix,levs ! ! t2m,q2m - real, 2 meter temperature and humidity im ! ! u10m,v10m - real, 10 meater u/v wind speed im ! ! zlvl - real, layer 1 height (m) im ! ! psurf - real, surface pressure (Pa) im ! ! hpbl - real, pbl height (m) im ! ! pwat - real, precipitable water im ! ! t1 - real, layer 1 temperature (K) im ! ! q1 - real, layer 1 specific humidity (kg/kg) im ! ! u1 - real, layer 1 zonal wind (m/s) im ! ! v1 - real, layer 1 merdional wind (m/s) im ! ! chh - real, thermal exchange coefficient im ! ! cmm - real, momentum exchange coefficient im ! ! dlwsfci - real, instantaneous sfc dnwd lw flux ( w/m**2 ) im ! ! ulwsfci - real, instantaneous sfc upwd lw flux ( w/m**2 ) im ! ! dswsfci - real, instantaneous sfc dnwd sw flux ( w/m**2 ) im ! ! uswsfci - real, instantaneous sfc upwd sw flux ( w/m**2 ) im ! ! dusfci - real, instantaneous u component of surface stress im ! ! dvsfci - real, instantaneous v component of surface stress im ! ! dtsfci - real, instantaneous sfc sensible heat flux im ! ! dqsfci - real, instantaneous sfc latent heat flux im ! ! gfluxi - real, instantaneous sfc ground heat flux im ! ! epi - real, instantaneous sfc potential evaporation im ! ! smcwlt2 - real, wilting point (volumetric) im ! ! smcref2 - real, soil moisture threshold (volumetric) im ! ! ! dusfci_cpl - real, sfc u-momentum flux at time step AOI cpl im ! ! dvsfci_cpl - real, sfc v-momentum flux at time step AOI cpl im ! ! dtsfci_cpl - real, sfc sensib heat flux at time step AOI cpl im ! ! dqsfci_cpl - real, sfc latent heat flux at time step AOI cpl im ! ! dlwsfci_cpl - real, sfc dnwd lw flux at time step AOI cpl im ! ! dswsfci_cpl - real, sfc dnwd sw flux at time step AOI cpl im ! ! dnirbmi_cpl - real, sfc nir beam dnwd sw flx rad at time step im ! ! dnirdfi_cpl - real, sfc nir diff dnwd sw flx rad at time step im ! ! dvisbmi_cpl - real, sfc uv+vis beam dnwd sw flx at time step im ! ! dvisdfi_cpl - real, sfc uv+vis diff dnwd sw flx at time step im ! ! nlwsfci_cpl - real, net sfc dnwd lw flux at time step AOI cpl im ! ! nswsfci_cpl - real, net sfc dnwd sw flux at time step AOI cpl im ! ! nnirbmi_cpl - real, net nir beam dnwd sw flx rad at time step im ! ! nnirdfi_cpl - real, net nir diff dnwd sw flx rad at time step im ! ! nvisbmi_cpl - real, net uv+vis beam dnwd sw flx at time step im ! ! nvisdfi_cpl - real, net uv+vis diff dnwd sw flx at time step im ! ! ocalnirbm_cpl- real, ocean alb nir beam (no ice) at time step im ! ! ocalnirdf_cpl- real, ocean alb nir diff (no ice) at time step im ! ! ocalvisbm_cpl- real, ocean alb vis beam (no ice) at time step im ! ! ocalvisdf_cpl- real, ocean alb vis diff (no ice) at time step im ! ! t2mi_cpl - real, T2m at time step AOI cpl im ! ! q2mi_cpl - real, Q2m at time step AOI cpl im ! ! u10mi_cpl - real, U10m at time step AOI cpl im ! ! v10mi_cpl - real, V10m at time step AOI cpl im ! ! tseai_cpl - real, sfc temp at time step AOI cpl im ! ! psurfi_cpl - real, sfc pressure at time step AOI cpl im ! ! tref - real, Reference Temperature im ! ! z_c - real, Sub-layer cooling thickness im ! ! c_0 - real, coefficient1 to calculate d(Tz)/d(Ts) im ! ! c_d - real, coefficient2 to calculate d(Tz)/d(Ts) im ! ! w_0 - real, coefficient3 to calculate d(Tz)/d(Ts) im ! ! w_d - real, coefficient4 to calculate d(Tz)/d(Ts) im ! ! rqtk - real, mass change due to moisture variation im ! ! dtdtr - real, temperature change due to radiative heating ! ! per time step (K) ix,levs! ! ! ! ==================== end of description ===================== ! subroutine gbphys & ! --- inputs: & ( im,ix,levs,lsoil,lsm,ntrac,ncld,ntoz,ntcw,ntke, & & nmtvr,nrcm,ko3,lonr,latr,jcap, & & num_p3d,num_p2d,npdf3d,ncnvcld3d, & & kdt,lat,me,pl_coeff,nlons,ncw,flgmin,crtrh,cdmbgwd, & & ccwf,dlqf,ctei_rm,clstp,cgwf,prslrd0,ral_ts,dtp,dtf,fhour, & & solhr,slag,sdec,cdec,sinlat,coslat,pgr,ugrs,vgrs, & & tgrs,qgrs,vvel,prsi,prsl,prslk,prsik,phii,phil, & & rann,prdout,poz,dpshc,fscav,fswtr,hprime,xlon,xlat, & & slope,shdmin,shdmax,snoalb,tg3,slmsk,vfrac, & & vtype,stype,uustar,oro,oro_uf,coszen,sfcdsw,sfcnsw, & & sfcnirbmd,sfcnirdfd,sfcvisbmd,sfcvisdfd, & & sfcnirbmu,sfcnirdfu,sfcvisbmu,sfcvisdfu, & & slimskin_cpl,ulwsfcin_cpl, & & dusfcin_cpl,dvsfcin_cpl,dtsfcin_cpl,dqsfcin_cpl, & & sfcdlw,tsflw,sfcemis,sfalb,swh,swhc,hlw,hlwc,hlwd,lsidea, & & ras,pre_rad,ldiag3d,lgocart,lssav,lssav_cpl, & & xkzm_m,xkzm_h,xkzm_s,psautco,prautco,evpco,wminco, & & pdfcld,shcnvcw,sup,redrag,hybedmf,dspheat, & & flipv,old_monin,cnvgwd,shal_cnv, & & imfshalcnv,imfdeepcnv,cal_pre, & & mom4ice,mstrat,trans_trac,nstf_name,moist_adj, & & thermodyn_id, sfcpress_id, gen_coord_hybrid,levr,adjtrc,nnp,& & cscnv,nctp,do_shoc,shocaftcnv,ntot3d,ntot2d, & ! --- input/outputs: & hice,fice,tisfc,tsea,tprcp,cv,cvb,cvt, & & srflag,snwdph,weasd,sncovr,zorl,canopy, & & ffmm,ffhh,f10m,srunoff,evbsa,evcwa,snohfa, & & transa,sbsnoa,snowca,soilm,tmpmin,tmpmax, & & dusfc,dvsfc,dtsfc,dqsfc,totprcp,gflux, & & dlwsfc,ulwsfc,suntim,runoff,ep,cldwrk, & & dugwd,dvgwd,psmean,cnvprcp,spfhmin,spfhmax,rain,rainc, & & dt3dt,dq3dt,du3dt,dv3dt,dqdt_v,cnvqc_v,acv,acvb,acvt, & & slc,smc,stc,upd_mf,dwn_mf,det_mf,phy_f3d,phy_f2d, & ! & slc,smc,stc,upd_mf,dwn_mf,det_mf,dkh,rnp,phy_f3d,phy_f2d, & & dusfc_cpl, dvsfc_cpl, dtsfc_cpl, dqsfc_cpl, & & dlwsfc_cpl,dswsfc_cpl,dnirbm_cpl,dnirdf_cpl, & & dvisbm_cpl,dvisdf_cpl,rain_cpl, nlwsfc_cpl,nswsfc_cpl, & & nnirbm_cpl,nnirdf_cpl,nvisbm_cpl,nvisdf_cpl,snow_cpl, & & xt,xs,xu,xv,xz,zm,xtts,xzts,d_conv,ifd,dt_cool,Qrain, & & phy_fctd, & ! --- outputs: & gt0,gq0,gu0,gv0,t2m,q2m,u10m,v10m, & & zlvl,psurf,hpbl,pwat,t1,q1,u1,v1, & & chh,cmm,dlwsfci,ulwsfci,dswsfci,uswsfci,dusfci,dvsfci, & & dtsfci,dqsfci,gfluxi,epi,smcwlt2,smcref2,wet1,sr, & & rqtk, & ! Stochastic physics perturnbation & dtdtr, & & dusfci_cpl, dvsfci_cpl, dtsfci_cpl, dqsfci_cpl, & & dlwsfci_cpl, dswsfci_cpl, & & dnirbmi_cpl, dnirdfi_cpl, dvisbmi_cpl, dvisdfi_cpl, & & nlwsfci_cpl, nswsfci_cpl, & & nnirbmi_cpl, nnirdfi_cpl, nvisbmi_cpl, nvisdfi_cpl, & & t2mi_cpl, q2mi_cpl, u10mi_cpl, v10mi_cpl, & & tseai_cpl, psurfi_cpl, & & tref, z_c, c_0, c_d, w_0, w_d & & ) ! use machine , only : kind_phys use physcons, only : con_cp, con_fvirt, con_g, con_rd, con_rv, & & con_hvap, con_hfus, con_rerth, con_pi &, rhc_max, dxmin, dxinv, pa2mb, rlapse use module_nst_water_prop, only: get_dtzm_2d use cs_conv, only : cs_convr implicit none ! ! --- some constant parameters: real(kind=kind_phys), parameter :: hocp = con_hvap/con_cp ! &, fhourpr = 0.0 &, qmin = 1.0e-10 &, p850 = 85000.0 &, epsq = 1.e-20 &, hsub = con_hvap+con_hfus &, czmin = 0.0001 ! cos(89.994) ! --- inputs: ! note: lgocart is the logical flag for in-line gocart; integer, intent(in) :: ix, im, levs, lsoil, lsm, ntrac, & & ncld, ntoz, ntcw, nmtvr, nrcm, ko3, & & lonr, latr, jcap, num_p3d, num_p2d, kdt, & & me, pl_coeff, lat, npdf3d, ncnvcld3d, & & thermodyn_id, sfcpress_id, levr, nnp, nctp,& & ntke, ntot3d, ntot2d integer, intent(in) :: nlons(im), ncw(2) integer, intent(in) :: nstf_name(5) integer, intent(in) :: imfshalcnv, imfdeepcnv logical, intent(in) :: ras, pre_rad, ldiag3d, flipv, & & old_monin, cnvgwd, & & redrag, hybedmf, dspheat, & & lssav, mom4ice, mstrat, & & trans_trac, moist_adj, cal_pre, cscnv, & & shal_cnv, gen_coord_hybrid, lgocart, & & lsidea, lssav_cpl, pdfcld, shcnvcw, & & do_shoc, shocaftcnv real(kind=kind_phys) :: adjtrc(ntrac) real(kind=kind_phys), dimension(im), intent(in) :: & & sinlat, coslat, pgr, dpshc, xlon, xlat, & & slope, shdmin, shdmax, snoalb, tg3, slmsk, vfrac, & & vtype, stype, uustar, oro, coszen, sfcnsw, sfcdsw, & & sfcnirbmu, sfcnirdfu, sfcvisbmu, sfcvisdfu, & & sfcnirbmd, sfcnirdfd, sfcvisbmd, sfcvisdfd, & & slimskin_cpl, dusfcin_cpl, dvsfcin_cpl, & & dtsfcin_cpl, dqsfcin_cpl, ulwsfcin_cpl, & & sfcdlw, tsflw, sfalb, sfcemis, oro_uf real(kind=kind_phys), dimension(ix,levs), intent(in) :: & & ugrs, vgrs, tgrs, vvel, prsl, prslk, phil, swh, swhc, hlw, hlwc !idea add by hmhj real(kind=kind_phys), intent(in) :: hlwd(ix,levs,6) real(kind=kind_phys), intent(inout) :: qgrs(ix,levs,ntrac) real(kind=kind_phys), dimension(ix,levs+1), intent(in) :: & & prsi, prsik, phii real(kind=kind_phys), intent(in) :: hprime(ix,nmtvr), & & prdout(ix,ko3,pl_coeff), rann(ix,nrcm), poz(ko3) real(kind=kind_phys), intent(in) :: dtp, dtf, fhour, solhr, & & slag, sdec, cdec, ctei_rm(2), clstp, & & ccwf(2), crtrh(3), flgmin(2), dlqf(2), cdmbgwd(2), & & xkzm_m, xkzm_h, xkzm_s, psautco(2), prautco(2), & & evpco, wminco(2), cgwf(2), prslrd0, sup, ral_ts ! --- input/output: real(kind=kind_phys), dimension(im), intent(inout) :: & & hice, fice, tisfc, tsea, tprcp, cv, cvb, cvt, & & srflag, snwdph, weasd, sncovr, zorl, canopy, ffmm, ffhh,& & f10m, srunoff, evbsa, evcwa, snohfa, transa, sbsnoa, & & snowca, soilm, tmpmin, tmpmax, dusfc, dvsfc, dtsfc, & & dqsfc, totprcp, gflux, dlwsfc, ulwsfc, suntim, runoff, ep,& & cldwrk, dugwd, dvgwd, psmean, cnvprcp,spfhmin,spfhmax, & & rain, rainc, acv, acvb, acvt real(kind=kind_phys), dimension(im), optional, intent(inout) :: & ! for A/O/I coupling & dusfc_cpl, dvsfc_cpl, dtsfc_cpl, dqsfc_cpl, & & dlwsfc_cpl,dswsfc_cpl,rain_cpl,snow_cpl, & & dnirbm_cpl,dnirdf_cpl,dvisbm_cpl,dvisdf_cpl, & & nlwsfc_cpl,nswsfc_cpl, & & nnirbm_cpl,nnirdf_cpl,nvisbm_cpl,nvisdf_cpl, & ! for nst & xt, xs, xu, xv, xz, zm, xtts, xzts, d_conv, ifd, dt_cool, & Qrain ! & Qrain, tref, z_c, c_0, c_d, w_0, w_d ! real(kind=kind_phys), dimension(ix,lsoil), intent(inout) :: & & smc, stc, slc real(kind=kind_phys), dimension(ix,levs), intent(inout) :: & & upd_mf, dwn_mf, det_mf, dqdt_v, cnvqc_v ! & upd_mf, dwn_mf, det_mf, dkh, rnp real(kind=kind_phys), intent(inout) :: & & phy_f3d(ix,levs,ntot3d), phy_f2d(ix,ntot2d), & & dt3dt(ix,levs,6), du3dt(ix,levs,4), dv3dt(ix,levs,4), & & dq3dt(ix,levs,5+pl_coeff) real(kind=kind_phys), intent(inout) :: & & phy_fctd(ix,nctp) & real(kind=kind_phys), dimension(ntrac-ncld+2) :: fscav, fswtr ! --- output: real(kind=kind_phys), dimension(im), intent(out) :: & & t2m, q2m, u10m, v10m, zlvl, psurf, hpbl, & & pwat, t1, q1, u1, v1, chh, cmm, & & dlwsfci, ulwsfci, dswsfci, uswsfci, & & dusfci, dvsfci, dtsfci, dqsfci, & & gfluxi, epi, smcwlt2, smcref2, wet1, sr real(kind=kind_phys), dimension(im), optional, intent(out) :: & & dusfci_cpl,dvsfci_cpl,dtsfci_cpl,dqsfci_cpl, & & dlwsfci_cpl,dswsfci_cpl, & & dnirbmi_cpl,dnirdfi_cpl,dvisbmi_cpl,dvisdfi_cpl, & & nlwsfci_cpl,nswsfci_cpl, & & nnirbmi_cpl,nnirdfi_cpl,nvisbmi_cpl,nvisdfi_cpl, & & t2mi_cpl,q2mi_cpl, & & u10mi_cpl,v10mi_cpl,tseai_cpl,psurfi_cpl, & ! & rqtk & tref, z_c, c_0, c_d, w_0, w_d, rqtk real(kind=kind_phys), dimension(ix,levs), intent(out) :: & & gt0, gu0, gv0 real(kind=kind_phys), dimension(ix,levs,ntrac), intent(out) :: & & gq0 ! --- local: ! real(kind=kind_phys) :: fscav(ntrac-ncld-1) ! real(kind=kind_phys),allocatable :: fscav(:), fswtr(:) real(kind=kind_phys), dimension(im) :: ccwfac, garea, & & dlength, xncw, cumabs, qmax, cice, zice, tice, & ! & gflx, rain, rainc, rainl, rain1, raincs, evapc, & & gflx, rainl, rain1, raincs, & & snowmt, cd, cdq, qss, dusfcg, dvsfcg, dusfc1, & & dvsfc1, dtsfc1, dqsfc1, rb, rhscnpy, drain, cld1d, & & evap, hflx, stress, t850, ep1d, gamt, gamq, & & sigmaf, oc, theta, gamma, sigma, & & elvmax, wind, work1, work2, runof, xmu, & ! & elvmax, wind, work1, work2, runof, xmu, oro_land, & & fm10, fh2, tsurf, tx1, tx2, ctei_r, flgmin_l, & & evbs, evcw, trans, sbsno, snowc, & & adjsfcdsw, adjsfcnsw, adjsfcdlw, adjsfculw,gabsbdlw, & & adjnirbmu, adjnirdfu, adjvisbmu, adjvisdfu, & & adjnirbmd, adjnirdfd, adjvisbmd, adjvisdfd, & & xcosz, tseal, snohf, dlqfac, work3, ctei_rml, cldf, & & domr, domzr, domip, doms, psautco_l, prautco_l real(kind=kind_phys), dimension(im) :: ocalnirbm_cpl, & & ocalnirdf_cpl,ocalvisbm_cpl,ocalvisdf_cpl ! & dswsfc, radsl, & ! & dlwsf1, ulwsf1, xcosz, tseal, snohf, dlqfac, & ! & domr, domzr, domip, doms ! real(kind=kind_phys), dimension(ix,levs) :: ud_mf, dd_mf, & ! & dt_mf, del real(kind=kind_phys), dimension(ix,levs) :: del, dtdtr real(kind=kind_phys), dimension(im,levs-1) :: dkt real(kind=kind_phys), dimension(im,levs) :: rhc, dtdt, & & dudt, dvdt, gwdcu, gwdcv, dtdtc, dmdt, & ! & diagn1, diagn2, cuhr, cumchr, & & qr_col, fc_ice, rainp, ud_mf, dd_mf, dt_mf, prnum ! & qr_col, fc_ice, rainp, ud_mf, dd_mf, dt_mf, shoc_cld, prnum real(kind=kind_phys), dimension(im,lsoil) :: smsoil, stsoil, & & ai, bi, cci, rhsmc, zsoil, slsoil real(kind=kind_phys) :: zsea1,zsea2 real(kind=kind_phys), dimension(im) :: dtzm real(kind=kind_phys) :: rhbbot, rhbtop, rhpbl, frain, f_rain, & & f_ice, qi, qw, qr, wc, tem, tem1, tem2, sume, sumr, sumq, & & dqdt(im,levs,ntrac), oa4(im,4), clx(im,4), albbm, albdf, & & xcosz_loc ! in clw, the first two varaibles are cloud water and ice. ! from third to ntrac are convective transportable tracers, ! third being the ozone, when ntrac=3 (valid only with ras) real(kind=kind_phys), allocatable :: clw(:,:,:), qpl(:,:),qpi(:,:) integer, dimension(im) :: kbot, ktop, kcnv, soiltyp, vegtype, & & kpbl, slopetyp, kinver, lmh, levshc, islmsk integer :: i, nvdiff, kk, ic, k, n, ipr, lv, k1, iter, levshcm, & & tracers, trc_shft, tottracer, num2, num3 & &, nshocm, nshoc, ntk logical, dimension(im) :: flag_iter, flag_guess, invrsn real(kind=kind_phys), dimension(im) :: dtsfc_cice, & & dqsfc_cice, dusfc_cice, dvsfc_cice, ulwsfc_cice, tisfc_cice, & & tsea_cice, hice_cice, fice_cice integer, dimension(im) :: islmsk_cice logical, dimension(im) :: flag_cice logical :: lprnt real(kind=kind_phys), allocatable :: cnvc(:,:),cnvw(:,:) real(kind=kind_phys) eng0, eng1, dtshoc ! ! for CS-convection ! real(kind=kind_phys), parameter :: wcbmax1=2.8, wcbmax2=1.4 real(kind=kind_phys), parameter :: wcbmax1=2.5, wcbmax2=1.5 ! real(kind=kind_phys), parameter :: wcbmax1=1.4, wcbmax2=1.4 real(kind=kind_phys) wcbmax(im) ! real(kind=kind_phys) tf, tcr, tcrf ! parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) parameter (tf=258.16, tcr=273.16, tcrf=1.0/(tcr-tf)) ! ! !===> ... begin here ! ! The option to compute tracer scavenging in GSM is disabled ! do i=1, ntrac-ncld-1 ! fscav(i) = 0. ! enddo ! --- ... set up check print point (for debugging) ! !************************************************************************* lprnt = .false. ! lprnt = me == 0 .and. kdt < 10 ! lprnt = kdt >= 19 ! ipr = 1 ! lprnt = kdt >= 19 ! if (me == 0 .and. kdt < 5) ! &write(0,*)' In gbphys:', im,ix,levs,lsoil,lsm, ! & ntrac,ncld,ntoz,ntcw, ! & nmtvr,nrcm,ko3,lonr,latr,jcap,num_p3d,num_p2d,npdf3d, ! & kdt,lat,me,pl_coeff,ncw,flgmin,crtrh,cdmbgwd ! &,' ccwf=',ccwf,' dlqf=',dlqf,' ras=',ras, ! & ' evpco=',evpco,' wminco=',wminco,' levr=',levr ! ipr = 1 ! lprnt = kdt .gt. 0 ! do i = 1, im ! work1(1) = xlon(i) * 57.29578 ! if (work1(1) >= 180.0) work1(1) = work1(1) - 360.0 ! work2(1) = xlat(i) * 57.29578 ! print *,' me=',me,' work1=',work1(1),' work2=',work2(1),' i=',i ! lprnt = kdt > 4320 ! lprnt = kdt > 0 .and. abs(work1(1)-8.4375) < 0.5 & ! & .and. abs(work2(1)+1.4) < 0.5 ! lprnt = kdt >= 14 .and. lat == 43 ! lprnt = kdt >= 256 .and. abs(xlon(i)*57.29578-136.07) < 0.08 & ! & .and. abs(xlat(i)*57.29578+1.1) < 0.101 ! lprnt = kdt >= 0 .and. abs(xlon(i)*57.29578-21.5) < 0.501 & ! & .and. abs(xlat(i)*57.29578+11.8) < 0.501 ! lprnt = kdt >= 0 .and. abs(xlon(i)*57.29578-8.4375) < 0.501 & ! & .and. abs(xlat(i)*57.29578+1.4) < 0.501 ! lprnt = kdt >= 40 .and. abs(xlon(i)*57.29578-288.75) < 1.501 & ! & .and. abs(xlat(i)*57.29578+12.38) < 1.501 ! lprnt = kdt >= 0 .and. abs(xlon(i)*57.29578-135.0) < 0.201 & ! & .and. abs(xlat(i)*57.29578-10.476) < 0.201 ! lprnt = kdt >= 0 .and. abs(xlon(i)*57.29578-110.3) < 0.201 & ! & .and. abs(xlat(i)*57.29578-2.0) < 0.201 ! if (kdt == 10) ! & print *,' i=',i,' xlon=',xlon(i)*57.29578, & ! & ' xlat=',xlat(i)*57.29578,' i=',i,' me=',me ! if (lprnt) then ! ipr = i ! exit ! endif ! enddo ! write(0,*)' In GBPHYS LSIDEA=',lsidea ! lprnt = .false. ! if(lprnt) then ! write(0,*)' im=',im,' ix=',ix,' levs=',levs,' lsoil=',lsoil, & ! & ' ntrac=',ntrac,' ntoz=',ntoz,' ntcw=',ntcw,' me=',me, & ! & ' xlat=',xlat(ipr),' kdt=',kdt,' slmsk=',slmsk(ipr), & ! & ' ntke=',ntke,' num_p3d=',num_p3d,' xlon=',xlon(ipr) ! & ' tsea=',tsea(ipr),' tref=',tref(ipr),' dt_cool=',dt_cool(ipr),& ! & ' dt_warm=',2.0*xt(ipr)/xz(ipr),' nrcm=',nrcm,' xlon=', ! & xlon(ipr), & ! & ' dt_warm=',dt_warm(ipr),' nrcm=',nrcm,' xlon=',xlon(ipr), & ! & ' sfalb=',sfalb(ipr),' kdt=',kdt ! write(0,*) ' pgr=',pgr(ipr),' kdt=',kdt,' ipr=',ipr ! write(0,*)' ipr=',ipr,' phy_f2d=',phy_f2d(ipr,1:num_p2d) ! write(0,*),' ugrs=',ugrs(ipr,:) ! write(0,*)' vgrs=',vgrs(ipr,:) ! write(0,*)' tgrs=',tgrs(ipr,:),' kdt=',kdt,' ipr=',ipr ! &, ' xlon=',xlon(ipr),' xlat=',xlat(ipr) ! write(0,*)' qgrs=',qgrs(ipr,:,1) ! write(0,*)' ozg=',qgrs(ipr,:,2) ! write(0,*)' tke=',qgrs(ipr,:,4) ! print *,' clw=',qgrs(ipr,:,3) ! &, ' xlon=',xlon(ipr),' xlat=',xlat(ipr) ! endif ! !************************************************************************* ! nvdiff = ntrac ! vertical diffusion of all tracers! ! ! --- ... figure out number of extra tracers ! tottracer = 0 ! no convective transport of tracers if (trans_trac) then if (ntcw > 0) then if (ntoz < ntcw) then trc_shft = ntcw + ncld - 1 else trc_shft = ntoz endif elseif (ntoz > 0) then trc_shft = ntoz else trc_shft = 1 endif tracers = ntrac - trc_shft tottracer = tracers if (ntoz > 0) tottracer = tottracer + 1 ! ozone is added separately endif if (ntke > 0) ntk = ntke - trc_shft + 3 ! if (lprnt) write(0,*)' trans_trac=',trans_trac,' tottracer=', & ! write(0,*)' trans_trac=',trans_trac,' tottracer=', & ! & tottracer,' trc_shft=',trc_shft,' kdt=',kdt ! &, ntrac-ncld+2,' clstp=',clstp,' kdt=',kdt ! &,' ntk=',ntk,' lat=',lat allocate ( clw(ix,levs,tottracer+2) ) if (do_shoc) then allocate (qpl(im,levs), qpi(im,levs)) endif if (.not. ras .or. .not. cscnv) then allocate ( cnvc(ix,levs), cnvw(ix,levs)) endif ! allocate (fscav(tottracer+3), fswtr(tottracer+3)) ! The option to compute tracer scavenging in GSM is disabled do i=1, tottracer+3 fscav(i) = 0. fswtr(i) = 0. enddo if (nnp == 1) then do n=1,ntrac if (abs(1.0-adjtrc(n)) > 1.0e-7) then do k=1,levs do i=1,im qgrs(i,k,n) = qgrs(i,k,n) * adjtrc(n) enddo enddo endif enddo endif ! call get_prs(im,ix,levs,ntrac,tgrs,qgrs, & & thermodyn_id, sfcpress_id, & & gen_coord_hybrid, & & prsi,prsik,prsl,prslk,phii,phil,del) ! & prsi,prsik,prsl,prslk,phii,phil,del,lprnt) ! ! if (lprnt) then ! write(0,*)' prsi=',prsi(ipr,:) ! write(0,*)' prsik=',prsik(ipr,:),' me=',me,' kdt=',kdt ! write(0,*)' prslk=',prslk(ipr,:),' me=',me,' kdt=',kdt ! write(0,*)' phil=',phil(ipr,:),' me=',me,' kdt=',kdt,' ipr=',ipr ! write(0,*)' prsl=',prsl(ipr,:),' kdt=',kdt ! print *,' del=',del(ipr,:) ! endif ! rhbbot = crtrh(1) rhpbl = crtrh(2) rhbtop = crtrh(3) ! ! --- ... frain=factor for centered difference scheme correction of rain amount. frain = dtf / dtp do i = 1, im sigmaf(i) = max( vfrac(i),0.01 ) ! sigmaf(i) = max( vfrac(i),0.3 ) if (lsm == 0) sigmaf(i) = 0.5 + vfrac(i) * 0.5 islmsk(i) = nint(slmsk(i)) if (islmsk(i) == 2) then soiltyp(i) = 9 vegtype(i) = 13 slopetyp(i) = 9 !! clu: qa(slopetyp) else soiltyp(i) = int( stype(i)+0.5 ) vegtype(i) = int( vtype(i)+0.5 ) slopetyp(i) = int( slope(i)+0.5 ) !! clu: slope -> slopetyp endif ! --- ... xw: transfer ice thickness & concentration from global to local variables zice(i) = hice(i) cice(i) = fice(i) tice(i) = tisfc(i) if (lssav_cpl) then islmsk_cice(i) = nint(slimskin_cpl(i)) flag_cice(i) = (islmsk_cice(i) == 4) ulwsfc_cice(i) = ulwsfcin_cpl(i) dusfc_cice(i) = dusfcin_cpl(i) dvsfc_cice(i) = dvsfcin_cpl(i) dtsfc_cice(i) = dtsfcin_cpl(i) dqsfc_cice(i) = dqsfcin_cpl(i) tisfc_cice(i) = tisfc(i) tsea_cice(i) = tsea(i) fice_cice(i) = fice(i) hice_cice(i) = hice(i) endif work1(i) = (log(coslat(i) / (nlons(i)*latr)) - dxmin) * dxinv work1(i) = max(0.0, min(1.0,work1(i))) work2(i) = 1.0 - work1(i) psurf(i) = pgr(i) work3(i) = prsik(i,1) / prslk(i,1) tem1 = con_rerth * (con_pi+con_pi)*coslat(i)/nlons(i) tem2 = con_rerth * con_pi / latr garea(i) = tem1 * tem2 dlength(i) = sqrt( tem1*tem1+tem2*tem2 ) cldf(i) = cgwf(1)*work1(i) + cgwf(2)*work2(i) wcbmax(i) = wcbmax1*work1(i) + wcbmax2*work2(i) enddo ! if (lprnt) write(0,*)' in gbphys work1&2=',work1(ipr),work2(ipr) ! &,' dxmin=',dxmin,' dxinv=',dxinv,' dx=', ! & log(coslat(ipr) / (nlons(ipr)*latr)) ! &,' coslat=',coslat(ipr),' nlons=',nlons(ipr),' latr=',latr ! --- ... transfer soil moisture and temperature from global to local variables do k = 1, lsoil do i = 1, im smsoil(i,k) = smc(i,k) stsoil(i,k) = stc(i,k) slsoil(i,k) = slc(i,k) !! clu: slc -> slsoil enddo enddo do k = 1, levs do i = 1, im dudt(i,k) = 0. dvdt(i,k) = 0. dtdt(i,k) = 0. dtdtc(i,k) = 0. enddo enddo do n = 1, ntrac do k = 1, levs do i = 1, im dqdt(i,k,n) = 0. enddo enddo enddo ! --- ... initialize dtdt with heating rate from dcyc2 ! if (lprnt) then ! do ipr=1,im ! write(0,*)' before dcyc2: im=',im,' lsoil=',lsoil,' levs=', & ! & levs & ! &, ' sde=',sdec,' cdec=',cdec,' tsea=',tsea(ipr),' ipr=',ipr & ! &, ' lat=',lat,' me=',me,' kdt=',kdt & ! &, ' sfcdlw=',sfcdlw(ipr),' sfcnsw=',sfcnsw(ipr) ! print *,' hlw=',hlw(ipr,:),' me=',me,' lat=',lat,xlon(ipr) ! print *,' swh=',swh(ipr,:),' me=',me,' lat=',lat,xlon(ipr) ! enddo ! endif ! --- ... adjust mean radiation fluxes and heating rates to fit for ! faster model time steps. ! sw: using cos of zenith angle as scaling factor ! lw: using surface air skin temperature as scaling factor if (pre_rad) then call dcyc2t3_pre_rad & ! --- inputs: & ( solhr,slag,sdec,cdec,sinlat,coslat, & & xlon,coszen,tsea,tgrs(1,1),tgrs(1,1), & & sfcdsw,sfcnsw,sfcdlw,swh,hlw, & & sfcnirbmu,sfcnirdfu,sfcvisbmu,sfcvisdfu, & & sfcnirbmd,sfcnirdfd,sfcvisbmd,sfcvisdfd, & & ix, im, levs, & ! --- input/output: & dtdt, & ! --- outputs: & adjsfcdsw,adjsfcnsw,adjsfcdlw,adjsfculw,xmu,xcosz, & & adjnirbmu,adjnirdfu,adjvisbmu,adjvisdfu, & & adjnirbmd,adjnirdfd,adjvisbmd,adjvisdfd & & ) else call dcyc2t3 & ! --- inputs: & ( solhr,slag,sdec,cdec,sinlat,coslat, & & xlon,coszen,tsea,tgrs(1,1),tsflw,sfcemis, & & sfcdsw,sfcnsw,sfcdlw,swh,swhc,hlw,hlwc, & & sfcnirbmu,sfcnirdfu,sfcvisbmu,sfcvisdfu, & & sfcnirbmd,sfcnirdfd,sfcvisbmd,sfcvisdfd, & & ix, im, levs, & ! --- input/output: & dtdt,dtdtc, & ! --- outputs: & adjsfcdsw,adjsfcnsw,adjsfcdlw,adjsfculw,xmu,xcosz, & & adjnirbmu,adjnirdfu,adjvisbmu,adjvisdfu, & & adjnirbmd,adjnirdfd,adjvisbmd,adjvisdfd & & ) ! ! save temp change due to radiation - need for sttp stochastic physics !--------------------------------------------------------------------- do k=1,levs do i=1,im dtdtr(i,k) = dtdtr(i,k) + dtdtc(i,k)*dtf enddo enddo endif ! if (lsidea) then !idea jw do k = 1, levs do i = 1, im ! dtdt(i,k) = hlwd(i,k,2) dtdt(i,k) = 0. enddo enddo endif ! --- convert lw fluxes for land/ocean/sea-ice models ! note: for sw: adjsfcdsw and adjsfcnsw are zenith angle adjusted downward/net fluxes. ! for lw: adjsfcdlw is (sfc temp adjusted) downward fluxe with no emiss effect. ! adjsfculw is (sfc temp adjusted) upward fluxe including emiss effect. ! one needs to be aware that that the absorbed downward lw flux (used by land/ocean ! models as downward flux) is not the same as adjsfcdlw but a value reduced by ! the factor of emissivity. however, the net effects are the same when seeing ! it either above the surface interface or below. ! ! - flux above the interface used by atmosphere model: ! down: adjsfcdlw; up: adjsfculw = sfcemis*sigma*T**4 + (1-sfcemis)*adjsfcdlw ! net = up - down = sfcemis * (sigma*T**4 - adjsfcdlw) ! - flux below the interface used by lnd/oc/ice models: ! down: sfcemis*adjsfcdlw; up: sfcemis*sigma*T**4 ! net = up - down = sfcemis * (sigma*T**4 - adjsfcdlw) ! --- ... define the downward lw flux absorbed by ground do i = 1, im gabsbdlw(i) = sfcemis(i) * adjsfcdlw(i) enddo ! if( lsidea ) then ! idea : moved temp adjust to idea_phys ! print *,' in gbphys: lsidea is true ' ! DTDT = 0. ! endif if (lssav) then ! --- ... accumulate/save output variables ! --- ... sunshine duration time is defined as the length of time (in mdl output ! interval) that solar radiation falling on a plane perpendicular to the ! direction of the sun >= 120 w/m2 do i = 1, im if ( xcosz(i) >= czmin ) then ! zenth angle > 89.994 deg tem1 = adjsfcdsw(i) / xcosz(i) if ( tem1 >= 120.0 ) then suntim(i) = suntim(i) + dtf endif endif enddo ! --- ... sfc lw fluxes used by atmospheric model are saved for output do i = 1, im dlwsfc(i) = dlwsfc(i) + adjsfcdlw(i)*dtf if (lssav_cpl) then if (flag_cice(i)) adjsfculw(i) = ulwsfc_cice(i) endif ulwsfc(i) = ulwsfc(i) + adjsfculw(i)*dtf psmean(i) = psmean(i) + pgr(i)*dtf ! mean surface pressure enddo if (ldiag3d) then if( lsidea ) then do k = 1, levs do i = 1, im dt3dt(i,k,1) = dt3dt(i,k,1) + hlwd(i,k,1)*dtf dt3dt(i,k,2) = dt3dt(i,k,2) + hlwd(i,k,2)*dtf dt3dt(i,k,3) = dt3dt(i,k,3) + hlwd(i,k,3)*dtf dt3dt(i,k,4) = dt3dt(i,k,4) + hlwd(i,k,4)*dtf dt3dt(i,k,5) = dt3dt(i,k,5) + hlwd(i,k,5)*dtf dt3dt(i,k,6) = dt3dt(i,k,6) + hlwd(i,k,6)*dtf enddo enddo else do k = 1, levs do i = 1, im dt3dt(i,k,1) = dt3dt(i,k,1) + hlw(i,k)*dtf dt3dt(i,k,2) = dt3dt(i,k,2) + swh(i,k)*dtf*xmu(i) enddo enddo endif endif endif ! end if_lssav_block do i = 1, im kcnv(i) = 0 kinver(i) = levs invrsn(i) = .false. tx1(i) = 0.0 tx2(i) = 10.0 ctei_r(i) = 10.0 enddo ! Only used for old shallow convection with mstrat=.true. if (((imfshalcnv == 0 .and. shal_cnv) .or. old_monin) & & .and. mstrat) then do i = 1, im ctei_rml(i) = ctei_rm(1)*work1(i) + ctei_rm(2)*work2(i) enddo do k = 1, levs/2 do i = 1, im if (prsi(i,1)-prsi(i,k+1) < 0.35*prsi(i,1) & & .and. (.not. invrsn(i))) then tem = (tgrs(i,k+1)-tgrs(i,k)) / (prsl(i,k)-prsl(i,k+1)) if ((tem > 0.00010 .and. tx1(i) < 0.0) .or. & (tem-abs(tx1(i)) > 0.0 .and. tx2(i) < 0.0)) then invrsn(i) = .true. if (qgrs(i,k,1) > qgrs(i,k+1,1)) then tem1 = tgrs(i,k+1) + hocp*max(qgrs(i,k+1,1),qmin) tem2 = tgrs(i,k) + hocp*max(qgrs(i,k,1),qmin) tem1 = tem1 / prslk(i,k+1) - tem2 / prslk(i,k) ! --- ... (cp/l)(deltathetae)/(deltatwater) > ctei_rm -> conditon for CTEI ctei_r(i) = (1.0/hocp)*tem1/(qgrs(i,k+1,1)-qgrs(i,k,1)& & + qgrs(i,k+1,ntcw)-qgrs(i,k,ntcw)) else ctei_r(i) = 10 endif if ( ctei_rml(i) > ctei_r(i) ) then kinver(i) = k else kinver(i) = levs endif endif tx2(i) = tx1(i) tx1(i) = tem endif enddo enddo endif ! --- ... check print ! ipr = 1 ! if (lprnt) then ! write(0,*)' before progtm: im=',im,' lsoil=',lsoil & ! &, ' nvdiff=',nvdiff,' adjsfcnsw=',adjsfcnsw(ipr) & ! &, ' adjsfcdlw=',adjsfcdlw(ipr),'adjsfculw=',adjsfculw(ipr)& ! &, ' sfcemis=',sfcemis(ipr),' tsea2=',tsea(ipr) & ! &, ' ipr=',ipr,' me=',me,' lat=',lat,' xlon=',xlon(ipr) & ! &, ' kdt=',kdt ! write(0,*)' dtdth=',dtdt(ipr,:),' kdt=',kdt ! endif ! if (lprnt) write(0,*)' phil=',phil(1,1:5) ! --- ... lu: initialize flag_guess, flag_iter, tsurf do i = 1, im tsurf(i) = tsea(i) flag_guess(i) = .false. flag_iter(i) = .true. drain(i) = 0.0 ep1d(i) = 0.0 runof(i) = 0.0 hflx(i) = 0.0 evap(i) = 0.0 evbs(i) = 0.0 evcw(i) = 0.0 trans(i) = 0.0 sbsno(i) = 0.0 snowc(i) = 0.0 snohf(i) = 0.0 zlvl(i) = phil(i,1) / con_g smcwlt2(i) = 0.0 smcref2(i) = 0.0 enddo ! --- ... lu: iter-loop over (sfc_diff,sfc_drv,sfc_ocean,sfc_sice) do iter = 1, 2 ! --- ... surface exchange coefficients ! ! if (lprnt) write(0,*)' tsea=',tsea(ipr),' tsurf=',tsurf(ipr),iter call sfc_diff(im,pgr,ugrs,vgrs,tgrs,qgrs,zlvl, & & tsea,zorl,cd,cdq,rb, & & prsl(1,1),work3,islmsk, & & stress,ffmm,ffhh, & & uustar,wind,phy_f2d(1,num_p2d),fm10,fh2, & & sigmaf,vegtype,shdmax, & & tsurf, flag_iter, redrag) ! if (lprnt) write(0,*)' cdq=',cdq(ipr),' iter=',iter & ! &, ' wind=',wind(ipr),'phy_f2d=',phy_f2d(ipr,num_p2d),' ugrs=' & ! &, ugrs(ipr,1),' vgrs=',vgrs(ipr,1) ! --- ... lu: update flag_guess do i = 1, im if (iter == 1 .and. wind(i) < 2.0) then flag_guess(i) = .true. endif enddo if ( nstf_name(1) > 0 ) then do i = 1, im if ( islmsk(i) == 0 ) then tem = (oro(i)-oro_uf(i)) * rlapse tseal(i) = tsea(i) + tem tsurf(i) = tsurf(i) + tem endif enddo ! if (lprnt) write(0,*)' tseaz1=',tsea(ipr),' tref=',tref(ipr) ! &, ' dt_cool=',dt_cool(ipr),' dt_warm=',2.0*(xt(ipr)/xz(ipr) ! &, ' kdt=',kdt ! &, ' tgrs=',tgrs(ipr,1),' prsl=',prsl(ipr,1) ! &, ' work3=',work3(ipr),' kdt=',kdt call sfc_nst & & ( im,lsoil,pgr,ugrs,vgrs,tgrs,qgrs,tref,cd,cdq, & & prsl(1,1),work3,islmsk,xlon,sinlat,stress, & & sfcemis,gabsbdlw,adjsfcnsw,tprcp,dtf,kdt,solhr,xcosz, & & phy_f2d(1,num_p2d),flag_iter,flag_guess,nstf_name, & & lprnt,ipr, & ! --- Input/output & tseal,tsurf,xt,xs,xu,xv,xz,zm,xtts,xzts,dt_cool, & & z_c,c_0,c_d,w_0,w_d,d_conv,ifd,qrain, & ! --- outputs: & qss, gflx, cmm, chh, evap, hflx, ep1d) ! if (lprnt) print *,' tseaz2=',tseal(ipr),' tref=',tref(ipr), ! & ' dt_cool=',dt_cool(ipr),' dt_warm=',2.0*xt(ipr)/xz(ipr), ! & ' kdt=',kdt do i = 1, im if ( islmsk(i) == 0 ) then tsurf(i) = tsurf(i) - (oro(i)-oro_uf(i)) * rlapse endif enddo ! --- ... run nsst model ... --- if ( nstf_name(1) > 1 ) then zsea1 = 0.001*real(nstf_name(4)) zsea2 = 0.001*real(nstf_name(5)) call get_dtzm_2d(xt,xz,dt_cool,z_c,real(islmsk(:)), & zsea1,zsea2,im,1,dtzm) do i = 1, im if ( islmsk(i) == 0 ) then tsea(i) = max(271.2,tref(i) + dtzm(i)) & -(oro(i)-oro_uf(i))*rlapse endif enddo endif ! if (lprnt) print *,' tseaz2=',tsea(ipr),' tref=',tref(ipr), & ! & ' dt_cool=',dt_cool(ipr),' dt_warm=',dt_warm(ipr),' kdt=',kdt else ! --- ... surface energy balance over ocean call sfc_ocean & ! --- inputs: & ( im,pgr,ugrs,vgrs,tgrs,qgrs,tsea,cd,cdq, & & prsl(1,1),work3,islmsk,phy_f2d(1,num_p2d),flag_iter, & ! --- outputs: & qss,cmm,chh,gflx,evap,hflx,ep1d & & ) endif ! if ( nstf_name(1) > 0 ) then ! if (lprnt) write(0,*)' sfalb=',sfalb(ipr),' ipr=',ipr & ! &, ' weasd=',weasd(ipr),' snwdph=',snwdph(ipr) & ! &, ' tprcp=',tprcp(ipr),' kdt=',kdt,' iter=',iter ! &,' tseabefland=',tsea(ipr) ! --- ... surface energy balance over land ! if (lsm == 1) then ! noah lsm call ! if (lprnt) write(0,*)' tsead=',tsea(ipr),' tsurf=',tsurf(ipr),iter call sfc_drv & ! --- inputs: & ( im,lsoil,pgr,ugrs,vgrs,tgrs,qgrs,soiltyp,vegtype,sigmaf, & & sfcemis,gabsbdlw,adjsfcdsw,adjsfcnsw,dtf,tg3,cd,cdq, & & prsl(1,1),work3,zlvl,islmsk,phy_f2d(1,num_p2d),slopetyp, & & shdmin,shdmax,snoalb,sfalb,flag_iter,flag_guess, & ! --- in/outs: & weasd,snwdph,tsea,tprcp,srflag,smsoil,stsoil,slsoil, & & canopy,trans,tsurf, & ! --- outputs: & sncovr,qss,gflx,drain,evap,hflx,ep1d,runof, & & cmm,chh,evbs,evcw,sbsno,snowc,soilm,snohf, & & smcwlt2,smcref2,zorl,wet1 & & ) else ! osu lsm call call sfc_land & ! --- inputs: & ( im,lsoil,pgr,ugrs,vgrs,tgrs,qgrs,smsoil,soiltyp, & & sigmaf,vegtype,sfcemis,adjsfcdlw,adjsfcnsw,dtf, & & tg3,cd,cdq,prsl(1,1),work3,islmsk, & ! & zorl,tg3,cd,cdq,prsl(1,1),work3,islmsk, & & phy_f2d(1,num_p2d),flag_iter,flag_guess, & ! --- input/outputs: & weasd,tsea,tprcp,srflag,stsoil,canopy,tsurf, & ! --- outputs: & qss,snowmt,gflx,zsoil,rhscnpy,rhsmc, & & ai,bi,cci,drain,evap,hflx,ep1d,cmm,chh, & & evbs,evcw,trans,sbsno,snowc,soilm, & & snohf,smcwlt2,smcref2 & & ) endif ! if (lprnt) write(0,*)' tseabeficemodel =',tsea(ipr),' me=',me & ! &, ' kdt=',kdt ! --- ... surface energy balance over seaice if (lssav_cpl) then do i = 1, im if (flag_cice(i)) then islmsk (i) = islmsk_cice(i) endif enddo endif call sfc_sice & ! --- inputs: & ( im,lsoil,pgr,ugrs,vgrs,tgrs,qgrs,dtf, & & sfcemis,gabsbdlw,adjsfcnsw,adjsfcdsw,srflag, & & cd,cdq,prsl(1,1),work3,islmsk,phy_f2d(1,num_p2d), & & flag_iter,mom4ice,lsm, lprnt,ipr, & ! & flag_iter,mom4ice,lsm, & ! --- input/outputs: & zice,cice,tice,weasd,tsea,tprcp,stsoil,ep1d, & ! --- outputs: & snwdph,qss,snowmt,gflx,cmm,chh,evap,hflx & & ) if (lssav_cpl) then do i = 1, im if (flag_cice(i)) then islmsk(i) = nint(slmsk(i)) endif enddo call sfc_cice & ! --- inputs: & ( im,ugrs,vgrs,tgrs,qgrs,cd,cdq,prsl(1,1),work3, & & islmsk_cice,phy_f2d(1,num_p2d),flag_iter, & & dqsfc_cice,dtsfc_cice, & ! --- outputs: & qss,cmm,chh,evap,hflx & & ) endif ! --- ... lu: update flag_iter and flag_guess do i = 1, im flag_iter(i) = .false. flag_guess(i) = .false. if(islmsk(i) == 1 .and. iter == 1) then if (wind(i) < 2.0) flag_iter(i) = .true. elseif (islmsk(i) == 0 .and. iter == 1 & & .and. nstf_name(1) > 0) then if (wind(i) < 2.0) flag_iter(i) = .true. endif enddo enddo ! end iter_loop do i = 1, im epi(i) = ep1d(i) dlwsfci(i) = adjsfcdlw(i) ulwsfci(i) = adjsfculw(i) uswsfci(i) = adjsfcdsw(i) - adjsfcnsw(i) dswsfci(i) = adjsfcdsw(i) gfluxi(i) = gflx(i) t1(i) = tgrs(i,1) q1(i) = qgrs(i,1,1) u1(i) = ugrs(i,1) v1(i) = vgrs(i,1) enddo if (lsm == 0) then ! osu lsm call do i = 1, im sncovr(i) = 0.0 if (weasd(i) > 0.0) sncovr(i) = 1.0 enddo endif ! --- ... update near surface fields call sfc_diag(im,pgr,ugrs,vgrs,tgrs,qgrs, & & tsea,qss,f10m,u10m,v10m,t2m,q2m,work3, & & evap,ffmm,ffhh,fm10,fh2) do i = 1, im phy_f2d(i,num_p2d) = 0.0 enddo if (lssav_cpl) then do i = 1, im dlwsfci_cpl(i) = adjsfcdlw(i) dswsfci_cpl(i) = adjsfcdsw(i) dlwsfc_cpl(i) = dlwsfc_cpl(i) + adjsfcdlw(i)*dtf dswsfc_cpl(i) = dswsfc_cpl(i) + adjsfcdsw(i)*dtf dnirbmi_cpl(i) = adjnirbmd(i) dnirdfi_cpl(i) = adjnirdfd(i) dvisbmi_cpl(i) = adjvisbmd(i) dvisdfi_cpl(i) = adjvisdfd(i) dnirbm_cpl(i) = dnirbm_cpl(i) + adjnirbmd(i)*dtf dnirdf_cpl(i) = dnirdf_cpl(i) + adjnirdfd(i)*dtf dvisbm_cpl(i) = dvisbm_cpl(i) + adjvisbmd(i)*dtf dvisdf_cpl(i) = dvisdf_cpl(i) + adjvisdfd(i)*dtf nlwsfci_cpl(i) = adjsfcdlw(i) - adjsfculw(i) nlwsfc_cpl(i) = nlwsfc_cpl(i) + nlwsfci_cpl(i)*dtf t2mi_cpl(i) = t2m(i) q2mi_cpl(i) = q2m(i) u10mi_cpl(i) = u10m(i) v10mi_cpl(i) = v10m(i) tseai_cpl(i) = tsea(i) psurfi_cpl(i) = pgr(i) enddo ! --- estimate mean albedo for ocean point without ice cover and apply ! them to net SW heat fluxes albdf = 0.06 do i = 1, im if (islmsk(i) /= 1) then ! not a land point ! --- compute open water albedo xcosz_loc = max( 0.0, min( 1.0, xcosz(i) )) ocalnirdf_cpl(i) = 0.06 ocalnirbm_cpl(i) = max(albdf, 0.026/(xcosz_loc**1.7+0.065) & & + 0.15 * (xcosz_loc-0.1) * (xcosz_loc-0.5) & & * (xcosz_loc-1.0)) ocalvisdf_cpl(i) = 0.06 ocalvisbm_cpl(i) = ocalnirbm_cpl(i) nnirbmi_cpl(i) = adjnirbmd(i)-adjnirbmd(i)*ocalnirbm_cpl(i) nnirdfi_cpl(i) = adjnirdfd(i)-adjnirdfd(i)*ocalnirdf_cpl(i) nvisbmi_cpl(i) = adjvisbmd(i)-adjvisbmd(i)*ocalvisbm_cpl(i) nvisdfi_cpl(i) = adjvisdfd(i)-adjvisdfd(i)*ocalvisdf_cpl(i) else nnirbmi_cpl(i) = adjnirbmd(i) - adjnirbmu(i) nnirdfi_cpl(i) = adjnirdfd(i) - adjnirdfu(i) nvisbmi_cpl(i) = adjvisbmd(i) - adjvisbmu(i) nvisdfi_cpl(i) = adjvisdfd(i) - adjvisdfu(i) endif nswsfci_cpl(i) = nnirbmi_cpl(i) + nnirdfi_cpl(i) & & + nvisbmi_cpl(i) + nvisdfi_cpl(i) nswsfc_cpl(i) = nswsfc_cpl(i) + nswsfci_cpl(i)*dtf nnirbm_cpl(i) = nnirbm_cpl(i) + nnirbmi_cpl(i)*dtf nnirdf_cpl(i) = nnirdf_cpl(i) + nnirdfi_cpl(i)*dtf nvisbm_cpl(i) = nvisbm_cpl(i) + nvisbmi_cpl(i)*dtf nvisdf_cpl(i) = nvisdf_cpl(i) + nvisdfi_cpl(i)*dtf enddo endif if (lssav) then do i = 1, im gflux(i) = gflux(i) + gflx(i) * dtf evbsa(i) = evbsa(i) + evbs(i) * dtf evcwa(i) = evcwa(i) + evcw(i) * dtf transa(i) = transa(i) + trans(i) * dtf sbsnoa(i) = sbsnoa(i) + sbsno(i) * dtf snowca(i) = snowca(i) + snowc(i) * dtf snohfa(i) = snohfa(i) + snohf(i) * dtf ep(i) = ep(i) + ep1d(i) * dtf tmpmax(i) = max(tmpmax(i),t2m(i)) tmpmin(i) = min(tmpmin(i),t2m(i)) spfhmax(i) = max(spfhmax(i),q2m(i)) spfhmin(i) = min(spfhmin(i),q2m(i)) enddo endif !!!!!!!!!!!!!!!!!Commented by Moorthi on July 18, 2012 !!!!!!!!!!!!!!!!!!! ! do i = 1, im ! --- ... compute coefficient of evaporation in evapc ! ! if (evapc(i) > 1.0e0) evapc(i) = 1.0e0 ! --- ... over snow cover or ice or sea, coef of evap =1.0e0 ! if (weasd(i) > 0.0 .or. slmsk(i) /= 1.0) evapc(i) = 1.0e0 ! enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! --- ... Boundary Layer and Free atmospheic turbulence parameterization ! if (lprnt) write(0,*)' tsea3=',tsea(ipr),' slmsk=',slmsk(ipr) & ! &, ' kdt=',kdt,' evap=',evap(ipr) ! if (lprnt) write(0,*)' dtdtb=',dtdt(ipr,1:10) ! do i = 1, im ! if (islmsk(i) == 0) then ! oro_land(i) = 0.0 ! else ! oro_land(i) = oro(i) ! endif ! enddo ! write(0,*)' before monin clstp=',clstp,' kdt=',kdt,' lat=',lat if (do_shoc) then call moninshoc(ix,im,levs,ntrac,ntcw,dvdt,dudt,dtdt,dqdt, & & ugrs,vgrs,tgrs,qgrs,phy_f3d(1,1,ntot3d-1), & ! tkh & prnum,ntke, & & prsik(1,1),rb,zorl,u10m,v10m,ffmm,ffhh, & & tsea,hflx,evap,stress,wind,kpbl, & & prsi,del,prsl,prslk,phii,phil,dtp, & & dusfc1,dvsfc1,dtsfc1,dqsfc1,dkt,hpbl, & & kinver, xkzm_m, xkzm_h, xkzm_s, lprnt, ipr,me) else if (hybedmf) then call moninedmf(ix,im,levs,nvdiff,ntcw,dvdt,dudt,dtdt,dqdt, & & ugrs,vgrs,tgrs,qgrs,swh,hlw,xmu, & & prsik(1,1),rb,zorl,u10m,v10m,ffmm,ffhh, & & tsea,qss,hflx,evap,stress,wind,kpbl, & & prsi,del,prsl,prslk,phii,phil,dtp,dspheat, & & dusfc1,dvsfc1,dtsfc1,dqsfc1,hpbl,gamt,gamq,dkt, & & kinver, xkzm_m, xkzm_h, xkzm_s, lprnt, ipr) elseif (.not. old_monin) then call moninq(ix,im,levs,nvdiff,ntcw,dvdt,dudt,dtdt,dqdt, & & ugrs,vgrs,tgrs,qgrs,swh,hlw,xmu, & & prsik(1,1),rb,ffmm,ffhh, & & tsea,qss,hflx,evap,stress,wind,kpbl, & & prsi,del,prsl,prslk,phii,phil,dtp,dspheat, & & dusfc1,dvsfc1,dtsfc1,dqsfc1,hpbl,gamt,gamq,dkt, & & kinver, xkzm_m, xkzm_h, xkzm_s, lprnt, ipr) else if (mstrat) then call moninp1(ix,im,levs,nvdiff,dvdt,dudt,dtdt,dqdt, & & ugrs,vgrs,tgrs,qgrs, & & prsik(1,1),rb,ffmm,ffhh,tsea,qss,hflx,evap,stress,wind, & & kpbl,prsi,del,prsl,prslk,phii,phil,dtp, & & dusfc1,dvsfc1,dtsfc1,dqsfc1,hpbl,gamt,gamq,dkt, & & kinver, xkzm_m, xkzm_h) ! & kinver, oro_land, ctei_r, ctei_rm, xkzm_m, xkzm_h) else call moninp(ix,im,levs,nvdiff,dvdt,dudt,dtdt,dqdt, & & ugrs,vgrs,tgrs,qgrs, & & prsik(1,1),rb,ffmm,ffhh,tsea,qss,hflx,evap,stress,wind, & & kpbl,prsi,del,prsl,phii,phil,dtp, & & dusfc1,dvsfc1,dtsfc1,dqsfc1,hpbl,gamt,gamq,dkt, & & xkzm_m,xkzm_h) endif endif ! end if_hybedmf endif ! end if_do_shoc if (lssav_cpl) then do i = 1, im if (flag_cice(i)) then cice(i) = fice_cice(i) tsea(i) = tsea_cice(i) dusfc1(i) = dusfc_cice(i) dvsfc1(i) = dvsfc_cice(i) dqsfc1(i) = dqsfc_cice(i) dtsfc1(i) = dtsfc_cice(i) endif enddo endif ! if (lprnt) then ! write(0,*) ' dusfc1=',dusfc1(ipr),' kdt=',kdt,' lat=',lat ! write(0,*)' dtsfc1=',dtsfc1(ipr) ! write(0,*)' dqsfc1=',dqsfc1(ipr) ! write(0,*)' dtdt=',dtdt(ipr,1:10) ! print *,' dudtm=',dudt(ipr,:) ! endif ! --- ... coupling insertion if (lssav_cpl) then do i=1, im dusfc_cpl(i) = dusfc_cpl(i) + dusfc1(i)*dtf dvsfc_cpl(i) = dvsfc_cpl(i) + dvsfc1(i)*dtf dtsfc_cpl(i) = dtsfc_cpl(i) + dtsfc1(i)*dtf dqsfc_cpl(i) = dqsfc_cpl(i) + dqsfc1(i)*dtf dusfci_cpl(i) = dusfc1(i) dvsfci_cpl(i) = dvsfc1(i) dtsfci_cpl(i) = dtsfc1(i) dqsfci_cpl(i) = dqsfc1(i) enddo endif !-------------------------------------------------------lssav if loop ---------- if (lssav) then do i = 1, im dusfc(i) = dusfc(i) + dusfc1(i)*dtf dvsfc(i) = dvsfc(i) + dvsfc1(i)*dtf dtsfc(i) = dtsfc(i) + dtsfc1(i)*dtf dqsfc(i) = dqsfc(i) + dqsfc1(i)*dtf dusfci(i) = dusfc1(i) dvsfci(i) = dvsfc1(i) dtsfci(i) = dtsfc1(i) dqsfci(i) = dqsfc1(i) enddo ! if (lprnt) then ! write(0,*)' dusfc=',dusfc(ipr),' dusfc1=',dusfc1(ipr),' dtf=', ! & dtf,' kdt=',kdt,' lat=',lat ! endif if (ldiag3d) then if (lsidea) then do k = 1, levs do i = 1, im dt3dt(i,k,3) = dt3dt(i,k,3) + dtdt(i,k)*dtf enddo enddo else do k = 1, levs do i = 1, im tem = dtdt(i,k) - (hlw(i,k)+swh(i,k)*xmu(i)) dt3dt(i,k,3) = dt3dt(i,k,3) + tem*dtf enddo enddo endif do k = 1, levs do i = 1, im du3dt(i,k,1) = du3dt(i,k,1) + dudt(i,k) * dtf du3dt(i,k,2) = du3dt(i,k,2) - dudt(i,k) * dtf dv3dt(i,k,1) = dv3dt(i,k,1) + dvdt(i,k) * dtf dv3dt(i,k,2) = dv3dt(i,k,2) - dvdt(i,k) * dtf enddo enddo ! update dqdt_v to include moisture tendency due to vertical diffusion ! if (lgocart) then ! do k = 1, levs ! do i = 1, im ! dqdt_v(i,k) = dqdt(i,k,1) * dtf ! enddo ! enddo ! endif do k = 1, levs do i = 1, im tem = dqdt(i,k,1) * dtf dq3dt(i,k,1) = dq3dt(i,k,1) + tem ! dqdt_v(i,k) = tem enddo enddo if (ntoz > 0) then do k = 1, levs do i = 1, im dq3dt(i,k,5) = dq3dt(i,k,5) + dqdt(i,k,ntoz) * dtf enddo enddo endif endif endif ! end if_lssav !-------------------------------------------------------lssav if loop ---------- ! ! Orographic gravity wave drag parameterization ! --------------------------------------------- if (nmtvr == 14) then ! current operational - as of 2014 do i = 1, im oc(i) = hprime(i,2) enddo do k = 1, 4 do i = 1, im oa4(i,k) = hprime(i,k+2) clx(i,k) = hprime(i,k+6) enddo enddo do i = 1, im theta(i) = hprime(i,11) gamma(i) = hprime(i,12) sigma(i) = hprime(i,13) elvmax(i) = hprime(i,14) enddo elseif (nmtvr == 10) then do i = 1, im oc(i) = hprime(i,2) enddo do k = 1, 4 do i = 1, im oa4(i,k) = hprime(i,k+2) clx(i,k) = hprime(i,k+6) enddo enddo elseif (nmtvr == 6) then do i = 1, im oc(i) = hprime(i,2) enddo do k = 1, 4 do i = 1, im oa4(i,k) = hprime(i,k+2) clx(i,k) = 0.0 enddo enddo else oc = 0 ; oa4 = 0 ; clx = 0 ; theta = 0 ; gamma = 0 ; sigma = 0 elvmax = 0 endif ! end if_nmtvr ! write(0,*)' before gwd clstp=',clstp,' kdt=',kdt,' lat=',lat call gwdps(im, ix, im, levs, dvdt, dudt, dtdt, & & ugrs, vgrs, tgrs, qgrs, & & kpbl, prsi, del, prsl, prslk, & & phii, phil, dtp, & & kdt, hprime(1,1), oc, oa4, clx, & & theta,sigma,gamma,elvmax,dusfcg, dvsfcg, & & con_g,con_cp,con_rd,con_rv, lonr, nmtvr, cdmbgwd, & & me, lprnt,ipr) ! if (lprnt) print *,' dudtg=',dudt(ipr,:) if (lssav) then do i = 1, im dugwd(i) = dugwd(i) + dusfcg(i)*dtf dvgwd(i) = dvgwd(i) + dvsfcg(i)*dtf enddo ! if (lprnt) print *,' dugwd=',dugwd(ipr),' dusfcg=',dusfcg(ipr) ! if (lprnt) print *,' dvgwd=',dvgwd(ipr),' dvsfcg=',dvsfcg(ipr) if (ldiag3d) then do k = 1, levs do i = 1, im du3dt(i,k,2) = du3dt(i,k,2) + dudt(i,k) * dtf dv3dt(i,k,2) = dv3dt(i,k,2) + dvdt(i,k) * dtf dt3dt(i,k,2) = dt3dt(i,k,2) + dtdt(i,k) * dtf enddo enddo endif endif if( .not. lsidea .and. ral_ts > 0.0) then ! call rayleigh_damp_mesopause(im, ix, im, levs, dvdt, dudt, dtdt, ! & ugrs, vgrs, dtp, con_cp, levr, prsl, prslrd0) ! else call rayleigh_damp(im, ix, im, levs, dvdt, dudt, dtdt, ugrs, & vgrs, dtp, con_cp, levr, pgr, prsl, & prslrd0, ral_ts) endif do k = 1, levs do i = 1, im gt0(i,k) = tgrs(i,k) + dtdt(i,k) * dtp gu0(i,k) = ugrs(i,k) + dudt(i,k) * dtp gv0(i,k) = vgrs(i,k) + dvdt(i,k) * dtp enddo enddo do n = 1, ntrac do k = 1, levs do i = 1, im gq0(i,k,n) = qgrs(i,k,n) + dqdt(i,k,n) * dtp enddo enddo enddo if( lsidea ) then ! idea convective adjustment call ideaca_up(prsi,gt0,ix,im,levs+1) endif ! --- ... check print ! if (me == 0) then ! sumq = 0.0 ! do k = 1, levs ! do i = 1, im ! sumq = sumq + (dqdt(i,k,1)+dqdt(i,k,ntcw)) * del(i,k) ! enddo ! enddo ! sume = 0.0 ! do i = 1, im ! sume = sume + dqsfc1(i) ! enddo ! sumq = sumq * 1000.0 / con_g ! sume = sume / con_hvap ! print *,' after mon: sumq=',sumq,' sume=',sume, ' kdt=',kdt ! endif ! --- ... ozone physics if (ntoz > 0 .and. ntrac >= ntoz) then call ozphys(ix,im,levs,ko3,dtp,gq0(1,1,ntoz),gq0(1,1,ntoz) & &, gt0, poz, prsl, prdout, pl_coeff, del, ldiag3d & &, dq3dt(1,1,6), me) endif ! --- ... to side-step the ozone physics ! if (ntrac >= 2) then ! do k = 1, levs ! gq0(k,ntoz) = qgrs(k,ntoz) ! enddo ! endif ! if (lprnt) then ! print *,' levs=',levs,' jcap=',jcap,' dtp',dtp & ! *, ' slmsk=',slmsk(ilon,ilat),' kdt=',kdt ! print *,' rann=',rann,' ncld=',ncld,' iq=',iq,' lat=',lat ! print *,' pgr=',pgr ! print *,' del=',del(ipr,:) ! print *,' prsl=',prsl(ipr,:) ! print *,' prslk=',prslk(ipr,:) ! print *,' rann=',rann(ipr,1) ! write(0,*)' gt0=',gt0(ipr,:) & ! &, ' kdt=',kdt,' xlon=',xlon(ipr),' xlat=',xlat(ipr) ! print *,' dtdt=',dtdt(ipr,:) ! print *,' gu0=',gu0(ipr,:) ! print *,' gv0=',gv0(ipr,:) ! write(0,*)' gq0=',(gq0(ipr,k,1),k=1,levs),' lat=',lat ! print *,' gq1=',(gq0(ipr,k,ntcw),k=1,levs) ! print *,' vvel=',vvel ! endif if (ldiag3d) then do k = 1, levs do i = 1, im dtdt(i,k) = gt0(i,k) ! dqdt(i,k,1) = gq0(i,k,1) dudt(i,k) = gu0(i,k) dvdt(i,k) = gv0(i,k) enddo enddo elseif (cnvgwd) then do k = 1, levs do i = 1, im dtdt(i,k) = gt0(i,k) enddo enddo endif ! end if_ldiag3d/cnvgwd if (ldiag3d .or. lgocart) then do k = 1, levs do i = 1, im dqdt(i,k,1) = gq0(i,k,1) enddo enddo endif ! end if_ldiag3d/lgocart call get_phi(im,ix,levs,ntrac,gt0,gq0, & & thermodyn_id, sfcpress_id, & & gen_coord_hybrid, & & prsi,prsik,prsl,prslk,phii,phil) ! if (lprnt) then ! print *,' phii2=',phii(ipr,:) ! print *,' phil2=',phil(ipr,:) ! endif do k = 1, levs do i = 1, im clw(i,k,1) = 0.0 clw(i,k,2) = -999.9 enddo enddo if (.not. ras .or. .not. cscnv) then do k = 1, levs do i = 1, im cnvc(i,k) = 0.0 cnvw(i,k) = 0.0 enddo enddo endif ! write(0,*)' before cnv clstp=',clstp,' kdt=',kdt,' lat=',lat ! --- ... for convective tracer transport (while using ras) if (ras .or. cscnv) then if (tottracer > 0) then if (ntoz > 0) then do k=1,levs do i=1,im clw(i,k,3) = gq0(i,k,ntoz) enddo enddo if (tracers > 0) then do n = 1, tracers do k=1,levs do i=1,im clw(i,k,3+n) = gq0(i,k,n+trc_shft) enddo enddo enddo endif else do n = 1, tracers do k=1,levs do i=1,im clw(i,k,2+n) = gq0(i,k,n+trc_shft) enddo enddo enddo endif endif endif ! end if_ras do i = 1, im ktop(i) = 1 kbot(i) = levs enddo ! --- ... calling condensation/precipitation processes ! -------------------------------------------- if (ntcw > 0) then do k = 1, levs do i = 1, im tem = rhbbot - (rhbbot-rhbtop) * (1.0-prslk(i,k)) tem = rhc_max * work1(i) + tem * work2(i) rhc(i,k) = max(0.0, min(1.0,tem)) enddo enddo if (num_p3d == 3) then ! brad ferrier's microphysics ! --- ... algorithm to separate different hydrometeor species do k = 1, levs do i = 1, im wc = gq0(i,k,ntcw) qi = 0. qr = 0. qw = 0. f_ice = max(0.0, min(1.0, phy_f3d(i,k,1))) f_rain = max(0.0, min(1.0, phy_f3d(i,k,2))) qi = f_ice*wc qw = wc-qi if (qw > 0.0) then qr = f_rain*qw qw = qw-qr endif ! if (f_ice >= 1.0) then ! qi = wc ! elseif (f_ice <= 0.0) then ! qw = wc ! else ! qi = f_ice*wc ! qw = wc-qi ! endif ! if (qw > 0.0 .and. f_rain > 0.0) then ! if (f_rain >= 1.0) then ! qr = qw ! qw = 0.0 ! else ! qr = f_rain*qw ! qw = qw-qr ! endif ! endif qr_col(i,k) = qr ! clw(i,k) = qi + qw clw(i,k,1) = qi clw(i,k,2) = qw ! --- ... array to track fraction of "cloud" in the form of ice ! if (qi+qw > epsq) then ! fc_ice(i,k) = qi / (qi+qw) ! else ! fc_ice(i,k) = 0.0 ! endif enddo enddo elseif (num_p3d == 4) then ! zhao-carr microphysics do i = 1, im psautco_l(i) = psautco(1)*work1(i) + psautco(2)*work2(i) prautco_l(i) = prautco(1)*work1(i) + prautco(2)*work2(i) enddo do k = 1, levs do i = 1, im clw(i,k,1) = gq0(i,k,ntcw) enddo enddo endif ! end if_num_p3d else ! if_ntcw do i = 1, im psautco_l(i) = psautco(1)*work1(i) + psautco(2)*work2(i) prautco_l(i) = prautco(1)*work1(i) + prautco(2)*work2(i) enddo do k = 1, levs do i = 1, im rhc(i,k) = 1.0 enddo enddo endif ! end if_ntcw ! ! Call SHOC if do_shoc is true and shocaftcnv is false ! if (do_shoc .and. .not. shocaftcnv) then if (num_p3d == 4) then do k=1,levs do i=1,im qpl(i,k) = 0.0 qpi(i,k) = 0.0 tem = gq0(i,k,ntcw) & & * max(0.0, MIN(1.0, (TCR-gt0(i,k))*TCRF)) clw(i,k,1) = tem ! ice clw(i,k,2) = gq0(i,k,ntcw) - tem ! water enddo enddo endif ! dtshoc = 60.0 ! dtshoc = 120.0 ! dtshoc = dtp ! nshocm = (dtp/dtshoc) + 0.001 ! dtshoc = dtp / nshocm ! do nshoc=1,nshocm ! if (lprnt) write(1000+me,*)' before shoc tke=',clw(ipr,:,ntk), ! &' kdt=',kdt,' lat=',lat,'xlon=',xlon(ipr),' xlat=',xlat(ipr) ! call shoc(ix, im, 1, levs, levs+1, dtshoc, me, lat, & call shoc(ix, im, 1, levs, levs+1, dtp, me, lat, & & prsl(1,1), phii(1,1), phil(1,1), & & gu0(1,1),gv0(1,1), vvel(1,1), gt0(1,1), gq0(1,1,1), & ! & clw(1,1,1), clw(1,1,2), qpi, qpl, rhc, shoc_cld(1,1)& ! & clw(1,1,1), clw(1,1,2), qpi, qpl, shoc_cld(1,1) & & clw(1,1,1), clw(1,1,2), qpi, qpl, & ! & sgs_cld(1:im,1:levs) & & phy_f3d(1,1,ntot3d-2), clw(1,1,ntk), hflx, evap, & & prnum, phy_f3d(1,1,ntot3d-1), phy_f3d(1,1,ntot3d), & & lprnt, ipr) ! do k=1,levs ! do i=1,im ! sgs_cld(i,k) = sgs_cld(i,k) + shoc_cld(i,k) ! enddo ! enddo ! if (lprnt) write(1000+me,*)' after shoc tke=',clw(1,:,ntk), ! &' kdt=',kdt ! enddo ! ! do k=1,levs ! write(1000+me,*)' maxcld=',maxval(sgs_cld(1:im,k)), ! write(1000+me,*)' maxtkh=',maxval(phy_f3d(1:im,k,ntot3d-1)), ! &' k=',k,' kdt=',kdt,' lat=',lat ! enddo ! write(0,*)' aft shoc gt0=',gt0(1,:),' lat=',lat ! write(0,*)' aft shoc gq0=',gq0(1,:,1),' lat=',lat ! write(0,*)' aft shoc gu0=',gu0(1,:),' lat=',lat ! endif ! if(do_shoc) ! --- ... calling convective parameterization ! if (.not. ras .and. .not. cscnv) then if (imfdeepcnv == 1) then ! no random cloud top call sascnvn(im,ix,levs,jcap,dtp,del,prsl,pgr,phil, & & clw,gq0,gt0,gu0,gv0,cld1d, & & rain1,kbot,ktop,kcnv,islmsk, & & vvel,ncld,ud_mf,dd_mf,dt_mf,cnvw,cnvc) elseif (imfdeepcnv == 2) then call mfdeepcnv(im,ix,levs,dtp,del,prsl,pgr,phil, & & clw,gq0,gt0,gu0,gv0,cld1d, & & rain1,kbot,ktop,kcnv,islmsk,garea, & & vvel,ncld,ud_mf,dd_mf,dt_mf,cnvw,cnvc) ! if (lprnt) print *,' rain1=',rain1(ipr) elseif (imfdeepcnv == 0) then ! random cloud top call sascnv(im,ix,levs,jcap,dtp,del,prsl,pgr,phil, & & clw,gq0,gt0,gu0,gv0,cld1d, & & rain1,kbot,ktop,kcnv,islmsk, & & vvel,rann,ncld,ud_mf,dd_mf,dt_mf,cnvw,cnvc) ! if (lprnt) print *,' rain1=',rain1(ipr),' rann=',rann(ipr,1) endif else ! ras ! if (lprnt) print *,' calling ras for kdt=',kdt,' me=',me & ! &, ' lprnt=',lprnt if(cscnv) then ! Chikira-Sugiyama !DD ! fscav(:) = 0.0 ! fswtr(:) = 0.0 call cs_convr( & !DD & ix ,im ,levs , tottracer+3 , & !DD & gt0 ,gq0 ,rain1 , clw , & !DD & phil ,phii , & !DD & prsl ,prsi , & !DD & dtp ,dtf , & !DD & ud_mf ,dd_mf ,dt_mf , & !DD & gu0 ,gv0 ,fscav, fswtr, & !DD ! & phy_fctd, me ) !DD & moorthi & phy_fctd, me, wcbmax ) !DD & moorthi ! & phy_fctd ) !DD do i = 1,im !DD rain1(i) = rain1(i) * (dtp*0.001) !DD enddo else if (ccwf(1) >= 0.0 .or. ccwf(2) >= 0 ) then do i=1,im ccwfac(i) = ccwf(1)*work1(i) + ccwf(2)*work2(i) dlqfac(i) = dlqf(1)*work1(i) + dlqf(2)*work2(i) lmh(i) = levs enddo else do i=1,im ccwfac(i) = -999.0 dlqfac(i) = 0.0 lmh(i) = levs enddo endif ! if (lprnt) write(0,*) ' calling ras for kdt=',kdt,' me=',me & ! &, ' lprnt=',lprnt,' ccwfac=',ccwfac(ipr) call rascnv(im, ix, levs, dtp, dtf, rann & &, gt0, gq0, gu0, gv0, clw, tottracer, fscav & &, prsi, prsl, prsik, prslk, phil, phii & &, kpbl, cd, rain1, kbot, ktop, kcnv & &, phy_f2d(1,num_p2d), flipv, pa2mb & &, me, garea, lmh, ccwfac, nrcm, rhc & &, ud_mf, dd_mf, dt_mf, dlqfac, lprnt, ipr, kdt) endif cld1d = 0 if (lgocart) then do k = 1, levs do i = 1, im upd_mf(i,k) = upd_mf(i,k) + ud_mf(i,k) * frain dwn_mf(i,k) = dwn_mf(i,k) + dd_mf(i,k) * frain det_mf(i,k) = det_mf(i,k) + dt_mf(i,k) * frain cnvqc_v(i,k) = cnvqc_v(i,k) + (clw(i,k,1)+clw(i,k,2)- & & gq0(i,k,ntcw)) * frain enddo enddo endif ! if (lgocart) ! --- ... update the tracers due to convective transport if (tottracer > 0) then if (ntoz > 0) then ! for ozone do k=1,levs do i=1,im gq0(i,k,ntoz) = clw(i,k,3) enddo enddo if (tracers > 0) then ! for other tracers do n = 1, tracers do k=1,levs do i=1,im gq0(i,k,n+trc_shft) = clw(i,k,3+n) enddo enddo enddo endif else do n = 1, tracers do k=1,levs do i=1,im gq0(i,k,n+trc_shft) = clw(i,k,2+n) enddo enddo enddo endif endif endif ! end if_not_ras ! do i = 1, im rainc(i) = frain * rain1(i) enddo ! if (lssav) then do i = 1, im cldwrk(i) = cldwrk(i) + cld1d(i) * dtf cnvprcp(i) = cnvprcp(i) + rainc(i) enddo if (ldiag3d) then do k = 1, levs do i = 1, im dt3dt(i,k,4) = dt3dt(i,k,4) + (gt0(i,k)-dtdt(i,k)) * frain dq3dt(i,k,2) = dq3dt(i,k,2) + (gq0(i,k,1)-dqdt(i,k,1)) & & * frain du3dt(i,k,3) = du3dt(i,k,3) + (gu0(i,k)-dudt(i,k)) * frain dv3dt(i,k,3) = dv3dt(i,k,3) + (gv0(i,k)-dvdt(i,k)) * frain upd_mf(i,k) = upd_mf(i,k) + ud_mf(i,k) * (con_g*frain) dwn_mf(i,k) = dwn_mf(i,k) + dd_mf(i,k) * (con_g*frain) det_mf(i,k) = det_mf(i,k) + dt_mf(i,k) * (con_g*frain) enddo enddo endif ! if (ldiag3d) endif ! end if_lssav ! ! update dqdt_v to include moisture tendency due to deep convection if (lgocart) then do k = 1, levs do i = 1, im dqdt_v(i,k) = (gq0(i,k,1)-dqdt(i,k,1)) * frain upd_mf(i,k) = upd_mf(i,k) + ud_mf(i,k) * frain dwn_mf(i,k) = dwn_mf(i,k) + dd_mf(i,k) * frain det_mf(i,k) = det_mf(i,k) + dt_mf(i,k) * frain cnvqc_v(i,k) = cnvqc_v(i,k) + (clw(i,k,1)+clw(i,k,2)) & *frain enddo enddo endif ! if (lgocart) ! if(npdf3d == 3 .and. num_p3d == 4) then num2 = num_p3d + 2 num3 = num2 + 1 do k = 1, levs do i = 1, im phy_f3d(i,k,num2) = cnvw(i,k) phy_f3d(i,k,num3) = cnvc(i,k) enddo enddo else if(npdf3d == 0 .and. ncnvcld3d == 1) then num2 = num_p3d + 1 do k = 1, levs do i = 1, im phy_f3d(i,k,num2) = cnvw(i,k) enddo enddo endif ! write(0,*)' cnvgwd moninshobl clstp=',clstp,' kdt=',kdt, ! & ' lat=',lat ! if (cnvgwd) then ! call convective gravity wave drag ! --- ... calculate maximum convective heating rate qmax [k/s] ! cuhr = temperature change due to deep convection do i = 1, im ! qmax(i) = 0. cumabs(i) = 0.0 work3(i) = 0.0 enddo do k = 1, levs do i = 1, im ! cuhr(i,k) = (gt0(i,k)-dtdt(i,k)) / dtf ! cuhr(i,k) = (gt0(i,k)-dtdt(i,k)) / dtp ! moorthi ! cumchr(i,k) = 0. ! gwdcu(i,k) = 0. ! gwdcv(i,k) = 0. ! diagn1(i,k) = 0. ! diagn2(i,k) = 0. if (k >= kbot(i) .and. k <= ktop(i)) then ! qmax(i) = max(qmax(i),cuhr(i,k)) ! cumabs(i) = cuhr(i,k) + cumabs(i) cumabs(i) = cumabs(i) + (gt0(i,k)-dtdt(i,k)) * del(i,k) work3(i) = work3(i) + del(i,k) endif enddo enddo do i=1,im if (work3(i) > 0.0) cumabs(i) = cumabs(i) / (dtp*work3(i)) enddo ! do i = 1, im ! do k = kbot(i), ktop(i) ! do k1 = kbot(i), k ! cumchr(i,k) = cuhr(i,k1) + cumchr(i,k) ! enddo ! cumchr(i,k) = cumchr(i,k) / cumabs(i) ! enddo ! enddo ! --- ... begin check print ****************************************** ! if (lprnt) then ! if (kbot(ipr) <= ktop(ipr)) then ! write(*,*) 'kbot <= ktop for (lat,lon) = ', & ! & xlon(ipr)*57.29578,xlat(ipr)*57.29578 ! write(*,*) 'kcnv kbot ktop qmax dlength ',kcnv(ipr), & ! & kbot(ipr),ktop(ipr),(86400.*qmax(ipr)),dlength(ipr) ! write(*,9000) kdt !9000 format(/,3x,'k',5x,'cuhr(k)',4x,'cumchr(k)',5x, & ! & 'at kdt = ',i4,/) ! do k = ktop(ipr), kbot(ipr),-1 ! write(*,9010) k,(86400.*cuhr(ipr,k)),(100.*cumchr(ipr,k)) !9010 format(2x,i2,2x,f8.2,5x,f6.0) ! enddo ! endif ! if (fhour >= fhourpr) then ! print *,' before gwdc in gbphys start print' ! write(*,*) 'fhour ix im levs = ',fhour,ix,im,levs ! print *,'dtp dtf = ',dtp,dtf ! write(*,9100) !9100 format(//,14x,'pressure levels',// & ! & ' ilev',7x,'prsi',8x,'prsl',8x,'delp',/) ! k = levs + 1 ! write(*,9110) k,(10.*prsi(ipr,k)) !9110 format(i4,2x,f10.3) ! do k = levs, 1, -1 ! write(*,9120) k,(10.*prsl(ipr,k)),(10.*del(ipr,k)) ! write(*,9110) k,(10.*prsi(ipr,k)) ! enddo !9120 format(i4,12x,2(2x,f10.3)) ! write(*,9130) !9130 format(//,10x,'before gwdc in gbphys',//,' ilev',6x, & ! & 'ugrs',9x,'gu0',8x,'vgrs',9x,'gv0',8x, & ! & 'tgrs',9x,'gt0',8x,'gt0b',8x,'dudt',8x,'dvdt',/) ! do k = levs, 1, -1 ! write(*,9140) k,ugrs(ipr,k),gu0(ipr,k), & ! & vgrs(ipr,k),gv0(ipr,k), & ! & tgrs(ipr,k),gt0(ipr,k),dtdt(ipr,k), & ! & dudt(ipr,k),dvdt(ipr,k) ! enddo !9140 format(i4,9(2x,f10.3)) ! print *,' before gwdc in gbphys end print' ! endif ! endif ! end if_lprnt ! --- ... end check print ******************************************** call gwdc(im, ix, im, levs, lat, ugrs, vgrs, tgrs, qgrs, & & prsl, prsi, del, cumabs, ktop, kbot, kcnv,cldf, & & con_g,con_cp,con_rd,con_fvirt, dlength, & & lprnt, ipr, fhour, & & gwdcu, gwdcv,dusfcg,dvsfcg) ! if (lprnt) then ! if (fhour >= fhourpr) then ! print *,' after gwdc in gbphys start print' ! write(*,9131) !9131 format(//,10x,'after gwdc in gbphys',//,' ilev',6x, & ! & 'ugrs',9x,'gu0',8x,'vgrs',9x,'gv0',8x, & ! & 'tgrs',9x,'gt0',8x,'gt0b',7x,'gwdcu',7x,'gwdcv',/) ! do k = levs, 1, -1 ! write(*,9141) k,ugrs(ipr,k),gu0(ipr,k), & ! & vgrs(ipr,k),gv0(ipr,k), & ! & tgrs(ipr,k),gt0(ipr,k),dtdt(ipr,k), & ! & gwdcu(ipr,k),gwdcv(ipr,k) ! enddo !9141 format(i4,9(2x,f10.3)) ! print *,' after gwdc in gbphys end print' ! endif ! endif ! --- ... write out cloud top stress and wind tendencies if (lssav) then do i = 1, im dugwd(i) = dugwd(i) + dusfcg(i)*dtf dvgwd(i) = dvgwd(i) + dvsfcg(i)*dtf enddo if (ldiag3d) then do k = 1, levs do i = 1, im du3dt(i,k,4) = du3dt(i,k,4) + gwdcu(i,k) * dtf dv3dt(i,k,4) = dv3dt(i,k,4) + gwdcv(i,k) * dtf ! du3dt(i,k,2) = du3dt(i,k,2) + diagn1(i,k) * dtf ! dv3dt(i,k,2) = dv3dt(i,k,2) + diagn2(i,k) * dtf enddo enddo endif endif ! end if_lssav ! --- ... update the wind components with gwdc tendencies do k = 1, levs do i = 1, im eng0 = 0.5*(gu0(i,k)*gu0(i,k)+gv0(i,k)*gv0(i,k)) gu0(i,k) = gu0(i,k) + gwdcu(i,k) * dtp gv0(i,k) = gv0(i,k) + gwdcv(i,k) * dtp eng1 = 0.5*(gu0(i,k)*gu0(i,k)+gv0(i,k)*gv0(i,k)) gt0(i,k) = gt0(i,k) + (eng0-eng1)/(dtp*con_cp) enddo enddo ! if (lprnt) then ! if (fhour >= fhourpr) then ! print *,' after tendency gwdc in gbphys start print' ! write(*,9132) !9132 format(//,10x,'after tendency gwdc in gbphys',//,' ilev',6x,& ! & 'ugrs',9x,'gu0',8x,'vgrs',9x,'gv0',8x, & ! & 'tgrs',9x,'gt0',8x,'gt0b',7x,'gwdcu',7x,'gwdcv',/) ! do k = levs, 1, -1 ! write(*,9142) k,ugrs(ipr,k),gu0(ipr,k),vgrs(ipr,k), & ! & gv0(ipr,k),tgrs(ipr,k),gt0(ipr,k),dtdt(ipr,k), & ! & gwdcu(ipr,k),gwdcv(ipr,k) ! enddo !9142 format(i4,9(2x,f10.3)) ! print *,' after tendency gwdc in gbphys end print' ! endif ! endif endif ! end if_cnvgwd (convective gravity wave drag) ! write(0,*)' after cnvgwd clstp=',clstp,' kdt=',kdt, ! & ' lat=',lat if (ldiag3d) then do k = 1, levs do i = 1, im dtdt(i,k) = gt0(i,k) ! dqdt(i,k,1) = gq0(i,k,1) enddo enddo endif if (ldiag3d .or. lgocart) then do k = 1, levs do i = 1, im dqdt(i,k,1) = gq0(i,k,1) enddo enddo endif ! write(0,*)' before do_shoc shal clstp=',clstp,' kdt=',kdt, ! & ' lat=',lat if (.not. do_shoc) then if (shal_cnv) then ! Shallow convection parameterization ! ----------------------------------- if (imfshalcnv == 1) then ! opr option now at 2014 !----------------------- call shalcnv(im,ix,levs,jcap,dtp,del,prsl,pgr,phil, & & clw,gq0,gt0,gu0,gv0, & & rain1,kbot,ktop,kcnv,islmsk, & & vvel,ncld,hpbl,hflx,evap,ud_mf,dt_mf, & & cnvw,cnvc) if (shcnvcw .and. num_p3d == 4 .and. npdf3d == 3 ) then do k = 1, levs do i = 1, im phy_f3d(i,k,num2) = cnvw(i,k) phy_f3d(i,k,num3) = cnvc(i,k) !??? phy_f3d(i,k,num2) = phy_f3d(i,k,num2) + cnvw(i,k) !??? phy_f3d(i,k,num3) = phy_f3d(i,k,num3) + cnvc(i,k) enddo enddo else if(npdf3d == 0 .and. ncnvcld3d == 1) then num2 = num_p3d + 1 do k = 1, levs do i = 1, im phy_f3d(i,k,num2) = cnvw(i,k) enddo enddo endif do i = 1, im raincs(i) = frain * rain1(i) rainc(i) = rainc(i) + raincs(i) enddo if (lssav) then do i = 1, im cnvprcp(i) = cnvprcp(i) + raincs(i) enddo endif elseif (imfshalcnv == 2) then call mfshalcnv(im,ix,levs,dtp,del,prsl,pgr,phil, & & clw,gq0,gt0,gu0,gv0, & & rain1,kbot,ktop,kcnv,islmsk,garea, & & vvel,ncld,hpbl,ud_mf,dt_mf,cnvw,cnvc) if (shcnvcw .and. num_p3d == 4 .and. npdf3d == 3 ) then do k = 1, levs do i = 1, im phy_f3d(i,k,num2) = cnvw(i,k) phy_f3d(i,k,num3) = cnvc(i,k) !??? phy_f3d(i,k,num2) = phy_f3d(i,k,num2) + cnvw(i,k) !??? phy_f3d(i,k,num3) = phy_f3d(i,k,num3) + cnvc(i,k) enddo enddo else if(npdf3d == 0 .and. ncnvcld3d == 1) then num2 = num_p3d + 1 do k = 1, levs do i = 1, im phy_f3d(i,k,num2) = cnvw(i,k) enddo enddo endif do i = 1, im raincs(i) = frain * rain1(i) rainc(i) = rainc(i) + raincs(i) enddo if (lssav) then do i = 1, im cnvprcp(i) = cnvprcp(i) + raincs(i) enddo endif elseif (imfshalcnv == 0) then ! modified Tiedtke Shallow convecton !----------------------------------- do i = 1, im levshc(i) = 0 enddo do k = 2, levs do i = 1, im if (prsi(i,1)-prsi(i,k) <= dpshc(i)) levshc(i) = k enddo enddo levshcm = 1 do i = 1, im levshcm = max(levshcm, levshc(i)) enddo ! if (lprnt) print *,' levshcm=',levshcm,' gt0sh=',gt0(ipr,:) ! &, ' lat=',lat if (mstrat) then ! As in CFSv2 call shalcv(im,ix,levshcm,dtp,del,prsi,prsl,prslk,kcnv, & & gq0,gt0,levshc,phil,kinver,ctei_r,ctei_rml & &, lprnt,ipr) else call shalcvt3(im,ix,levshcm,dtp,del,prsi,prsl,prslk, & & kcnv,gq0,gt0) endif ! if (lprnt) print *,' levshcm=',levshcm,' gt0sha=',gt0(ipr,:) endif ! end if_imfshalcnv endif ! end if_shal_cnv if (lssav) then ! update dqdt_v to include moisture tendency due to shallow convection if (lgocart) then do k = 1, levs do i = 1, im tem = (gq0(i,k,1)-dqdt(i,k,1)) * frain dqdt_v(i,k) = dqdt_v(i,k) + tem enddo enddo endif if (ldiag3d) then do k = 1, levs do i = 1, im dt3dt(i,k,5) = dt3dt(i,k,5) + (gt0(i,k)-dtdt(i,k)) & * frain dq3dt(i,k,3) = dq3dt(i,k,3) + (gq0(i,k,1)-dqdt(i,k,1)) & & * frain dtdt(i,k) = gt0(i,k) dqdt(i,k,1) = gq0(i,k,1) enddo enddo endif endif ! end if_lssav ! do k = 1, levs do i = 1, im if (clw(i,k,2) <= -999.0) clw(i,k,2) = 0.0 enddo enddo ! if (lprnt) then ! write(0,*)' prsl=',prsl(ipr,:) ! write(0,*) ' del=',del(ipr,:) ! write(0,*) ' befshgt0=',gt0(ipr,:) ! write(0,*) ' befshgq0=',gq0(ipr,:,1) ! endif elseif (shocaftcnv) then ! if do_shoc is true and shocaftcnv is true call shoc if (clw(1,1,2) < -999.0) then ! if clw is not partitioned to ice and water do k=1,levs do i=1,im tem = gq0(i,k,ntcw) & & * max(0.0, MIN(1.0, (TCR-gt0(i,k))*TCRF)) clw(i,k,1) = tem ! ice clw(i,k,2) = gq0(i,k,ntcw) - tem ! water enddo enddo endif do k=1,levs do i=1,im qpl(i,k) = 0.0 qpi(i,k) = 0.0 enddo enddo ! dtshoc = 60.0 ! nshocm = (dtp/dtshoc) + 0.001 ! dtshoc = dtp / nshocm ! do nshoc=1,nshocm ! call shoc(im, 1, levs, levs+1, dtp, me, lat, & !! call shoc(im, 1, levs, levs+1, dtshoc, me, lat, & ! & prsl(1:im,:), phii (1:im,:), phil(1:im,:),& ! & gu0(1:im,:),gv0(1:im,:), vvel(1:im,:), gt0(1:im,:), & ! & gq0(1:im,:,1), & ! & clw(1:im,:,1), clw(1:im,:,2), qpi, qpl, sgs_cld(1:im,:)& ! &, gq0(1:im,:,ntke), & ! & phy_f3d(1:im,:,ntot3d-1), phy_f3d(1:im,:,ntot3d), & ! & lprnt, ipr, & ! & con_cp, con_g, con_hvap, con_hfus, con_hvap+con_hfus, & ! & con_rv, con_rd, con_pi, con_fvirt) ! call shoc(ix, im, 1, levs, levs+1, dtshoc, me, lat, & call shoc(ix, im, 1, levs, levs+1, dtp, me, lat, & & prsl(1,1), phii(1,1), phil(1,1), & & gu0(1,1),gv0(1,1), vvel(1,1), gt0(1,1), gq0(1,1,1), & & clw(1,1,1), clw(1,1,2), qpi, qpl, & & phy_f3d(1,1,ntot3d-2), gq0(1,1,ntke),hflx,evap, & & prnum, phy_f3d(1,1,ntot3d-1), phy_f3d(1,1,ntot3d), & & lprnt, ipr) ! ! do k=1,levs ! write(1000+me,*)' maxtkh=',maxval(phy_f3d(1:im,k,ntot3d-1)), ! &' k=',k,' kdt=',kdt,' lat=',lat ! enddo ! write(0,*)' aft shoc gt0=',gt0(1,:),' lat=',lat ! write(0,*)' aft shoc gq0=',gq0(1,:,1),' lat=',lat ! write(0,*)' aft shoc gu0=',gu0(1,:),' lat=',lat ! endif ! if( .not. do_shoc) ! ! if (lprnt) then ! write(0,*)' prsl=',prsl(ipr,:) ! write(0,*) ' del=',del(ipr,:) ! write(0,*) ' aftshgt0=',gt0(ipr,:) ! write(0,*) ' aftshgq0=',gq0(ipr,:,1) ! endif if (ntcw > 0) then ! microphysics if (num_p3d == 3) then ! call brad ferrier's microphysics do k = 1, levs do i = 1, im ! qi = clw(i,k)*fc_ice(i,k) ! qw = clw(i,k) - qi qi = clw(i,k,1) qw = clw(i,k,2) ! --- ... algorithm to combine different hydrometeor species ! gq0(i,k,ntcw) = max(epsq, qi+qw+qr_col(i,k)) gq0(i,k,ntcw) = qi + qw + qr_col(i,k) if (qi <= epsq) then phy_f3d(i,k,1) = 0. else phy_f3d(i,k,1) = qi/gq0(i,k,ntcw) endif if (qr_col(i,k) <= epsq) then phy_f3d(i,k,2) = 0. else phy_f3d(i,k,2) = qr_col(i,k) / (qw+qr_col(i,k)) endif enddo enddo elseif (num_p3d == 4) then ! if_num_p3d do k = 1, levs do i = 1, im gq0(i,k,ntcw) = clw(i,k,1) + clw(i,k,2) enddo enddo endif ! end if_num_p3d else ! if_ntcw do k = 1, levs do i = 1, im clw(i,k,1) = clw(i,k,1) + clw(i,k,2) enddo enddo endif ! end if_ntcw ! write(0,*)' bef cnvc90 clstp=',clstp,' kdt=',kdt,' lat=',lat call cnvc90(clstp, im, ix, rainc, kbot, ktop, levs, prsi, & & acv, acvb, acvt, cv, cvb, cvt) if (moist_adj) then ! moist convective adjustment ! --------------------------- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! if (me == 0) then ! sumq = 0.0 ! DO K=1,LEVS ! do i=1,im ! sumq = sumq - (gq0(i,k,1)+gq0(i,k,ntcw)) * del(i,k) ! enddo ! enddo ! endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! To call moist convective adjustment ! ! if (lprnt) then ! print *,' prsl=',prsl(ipr,:) ! print *,' del=',del(ipr,:) ! print *,' gt0b=',gt0(ipr,:) ! print *,' gq0b=',gq0(ipr,:,1) ! endif call mstcnv(im,ix,levs,dtp,gt0,gq0,prsl,del,prslk,rain1 &, gq0(1,1,ntcw), rhc, lprnt,ipr) ! if (lprnt) then ! print *,' rain1=',rain1(ipr),' rainc=',rainc(ipr) ! print *,' gt0a=',gt0(ipr,:) ! print *,' gq0a=',gq0(ipr,:,1) ! endif do i=1,im rainc(i) = rainc(i) + frain * rain1(i) enddo if(lssav) then do i=1,im cnvprcp(i) = cnvprcp(i) + rain1(i) * frain enddo ! update dqdt_v to include moisture tendency due to surface processes ! dqdt_v : instaneous moisture tendency (kg/kg/sec) ! if (lgocart) then ! do k=1,levs ! do i=1,im ! tem = (gq0(i,k,1)-dqdt(i,k,1)) * frain ! dqdt_v(i,k) = dqdt_v(i,k) + tem ! dqdt_v(i,k) = dqdt_v(i,k) / dtf ! enddo ! enddo ! endif if (ldiag3d) then do k=1,levs do i=1,im dt3dt(i,k,4) = dt3dt(i,k,4) + (gt0(i,k)-dtdt(i,k)) & * frain dq3dt(i,k,2) = dq3dt(i,k,2) + (gq0(i,k,1)-dqdt(i,k,1)) & * frain enddo enddo endif endif ! if (ldiag3d) then do k=1,levs do i=1,im dtdt(i,k) = gt0(i,k) dqdt(i,k,1) = gq0(i,k,1) enddo enddo endif ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! if (me == 0) then ! DO K=1,LEVS ! do i=1,im ! sumq = sumq + (gq0(i,k,1)+gq0(i,k,ntcw)) * del(i,k) ! enddo ! enddo ! sumr = 0.0 ! do i=1,im ! sumr = sumr + rain1(i) ! enddo ! sumq = sumq * 1000.0 / grav ! sumr = sumr *1000 ! print *,' after MCN: sumq=',sumq,' sumr=',sumr, ' kdt=',kdt ! endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! endif ! moist convective adjustment over ! dqdt_v : instaneous moisture tendency (kg/kg/sec) if (lgocart) then do k=1,levs do i=1,im dqdt_v(i,k) = dqdt_v(i,k) / dtf enddo enddo endif ! ! grid-scale condensation/precipitations and microphysics parameterization ! ------------------------------------------------------------------------ if (ncld == 0) then ! no cloud microphysics call lrgscl(ix,im,levs,dtp,gt0,gq0,prsl,del,prslk,rain1,clw) elseif (ncld == 1) then ! microphysics with single cloud condensate if (num_p3d == 3) then ! brad ferrier's microphysics ! --------------------------- do i = 1, im xncw(i) = ncw(1) * work1(i) + ncw(2) * work2(i) flgmin_l(i) = flgmin(1)* work1(i) + flgmin(2) * work2(i) enddo if (kdt == 1 .and. abs(xlon(1)) < 0.0001) then write(0,*)' xncw=',xncw(1),' rhc=',rhc(1,1),' work1=' & &, work1(1),' work2=',work2(1),' flgmin=',flgmin_l(1) & &, ' lon=',xlon(1) * 57.29578,' lat=',lat,' me=',me endif call gsmdrive(im, ix, levs, dtp, con_g, con_hvap, hsub, con_cp& &, me, lprnt, ipr & &, prsl, del, rhc, xncw, flgmin_l & &, gt0, gq0(1,1,1), gq0(1,1,ntcw) & &, phy_f3d(1,1,1), phy_f3d(1,1,2) & &, phy_f3d(1,1,3), rain1, sr) elseif (num_p3d == 4) then ! call zhao/carr/sundqvist microphysics if (npdf3d /= 3) then ! without pdf clouds ! if (lprnt) then ! write(0,*)' prsl=',prsl(ipr,:) ! write(0,*) ' del=',del(ipr,:) ! write(0,*) ' beflsgt0=',gt0(ipr,:),' kdt=',kdt ! write(0,*) ' beflsgq0=',gq0(ipr,:,1),' kdt=',kdt ! endif ! ------------------ if (do_shoc) then call precpd_shoc(im, ix, levs, dtp, del, prsl, & & gq0(1,1,1), gq0(1,1,ntcw), gt0, rain1, sr, & & rainp, rhc, psautco_l, prautco_l, evpco, & & wminco, phy_f3d(1,1,ntot3d-2), lprnt, ipr) ! & wminco, sgs_cld(1:im,1:levs), lprnt, ipr) ! & wminco, shoc_cld, lprnt, ipr) else call gscond(im, ix, levs, dtp, dtf, prsl, pgr, & & gq0(1,1,1), gq0(1,1,ntcw), gt0, & & phy_f3d(1,1,1), phy_f3d(1,1,2), phy_f2d(1,1), & & phy_f3d(1,1,3), phy_f3d(1,1,4), phy_f2d(1,2), & & rhc,lprnt, ipr) call precpd(im, ix, levs, dtp, del, prsl, & & gq0(1,1,1), gq0(1,1,ntcw), gt0, rain1, sr, & & rainp, rhc, psautco_l, prautco_l, evpco, & & wminco, lprnt, ipr) endif ! if (lprnt) then ! write(0,*)' prsl=',prsl(ipr,:) ! write(0,*) ' del=',del(ipr,:) ! write(0,*) ' aftlsgt0=',gt0(ipr,:),' kdt=',kdt ! write(0,*) ' aftlsgq0=',gq0(ipr,:,1),' kdt=',kdt ! endif else ! with pdf clouds ! --------------- call gscondp(im, ix, levs, dtp, dtf, prsl, pgr, & & gq0(1,1,1), gq0(1,1,ntcw), gt0, & & phy_f3d(1,1,1), phy_f3d(1,1,2), phy_f2d(1,1), & & phy_f3d(1,1,3), phy_f3d(1,1,4), phy_f2d(1,2), & & rhc,phy_f3d(1,1,num_p3d+1),sup,lprnt, & & ipr,kdt) call precpdp(im, ix, levs, dtp, del, prsl, pgr, & & gq0(1,1,1), gq0(1,1,ntcw), gt0, rain1,sr, & & rainp, rhc, phy_f3d(1,1,num_p3d+1), & & psautco_l, prautco_l, evpco, wminco, & & lprnt, ipr) endif ! end of grid-scale precip/microphysics options endif ! end if_num_p3d endif ! end if_ncld ! if (lprnt) write(0,*) ' rain1=',rain1(ipr),' rainc=',rainc(ipr) do i = 1, im rainl(i) = frain * rain1(i) rain(i) = rainc(i) + rainl(i) enddo if (cal_pre) then ! hchuang: add dominant precipitation type algorithm call calpreciptype(kdt,nrcm,im,ix,levs,levs+1,rann, & xlat,xlon,gt0,gq0,prsl,prsi,rain, & phii,num_p3d,tsea,sr,phy_f3d(1,1,3), ! input & domr,domzr,domip,doms) ! output ! ! if (lprnt) print*,'debug calpreciptype: DOMR,DOMZR,DOMIP,DOMS ' ! &,DOMR(ipr),DOMZR(ipr),DOMIP(ipr),DOMS(ipr) ! do i=1,im ! if (abs(xlon(i)*57.29578-114.0) .lt. 0.2 .and. ! & abs(xlat(i)*57.29578-40.0) .lt. 0.2) ! & print*,'debug calpreciptype: DOMR,DOMZR,DOMIP,DOMS ', ! & DOMR(i),DOMZR(i),DOMIP(i),DOMS(i) ! end do ! HCHUANG: use new precipitation type to decide snow flag for LSM snow accumulation do i=1,im if(doms(i) >0.0 .or. domip(i)>0.0)then srflag(i) = 1. else srflag(i) = 0. end if enddo endif if (lssav) then do i = 1, im totprcp(i) = totprcp(i) + rain(i) enddo if (ldiag3d) then do k = 1, levs do i = 1, im dt3dt(i,k,6) = dt3dt(i,k,6) + (gt0(i,k)-dtdt(i,k)) * frain dq3dt(i,k,4) = dq3dt(i,k,4) + (gq0(i,k,1)-dqdt(i,k,1)) & & * frain enddo enddo endif endif ! if (lprnt) write(0,*) ' bef estimate t850 for' ! --- ... estimate t850 for rain-snow decision do i = 1, im t850(i) = gt0(i,1) enddo do k = 1, levs-1 do i = 1, im if (prsl(i,k) > p850 .and. prsl(i,k+1) <= p850) then t850(i) = gt0(i,k) - (prsl(i,k)-p850) & & / (prsl(i,k)-prsl(i,k+1)) * (gt0(i,k)-gt0(i,k+1)) endif enddo enddo ! --- ... lu: snow-rain detection is performed in land/sice module if (cal_pre) then ! hchuang: new precip type algorithm defines srflag do i = 1, im tprcp(i) = rain(i) ! clu: rain -> tprcp enddo else do i = 1, im tprcp(i) = rain(i) ! clu: rain -> tprcp srflag(i) = 0. ! clu: default srflag as 'rain' (i.e. 0) if (t850(i) <= 273.16) then srflag(i) = 1. ! clu: set srflag to 'snow' (i.e. 1) endif enddo endif ! --- ... coupling insertion if (lssav_cpl) then do i = 1, im if (t850(i) > 273.16) then rain_cpl(i) = rain_cpl(i) + rain(i) else snow_cpl(i) = snow_cpl(i) + rain(i) endif enddo endif ! --- ... end coupling insertion if (lsm == 0) then ! for OSU land model ! --- ... wei: when calling osu lsm update soil moisture and canopy water ! after precipitation computaion do i = 1, im if (t850(i) <= 273.16 .and. islmsk(i) /= 0) then weasd(i) = weasd(i) + 1.e3*rain(i) tprcp(i) = 0. endif enddo call progt2(im,lsoil,rhscnpy,rhsmc,ai,bi,cci,smsoil, & & islmsk,canopy,tprcp,runof,snowmt, & & zsoil,soiltyp,sigmaf,dtf,me) ! --- ... wei: let soil liquid water equal to soil total water do k = 1, lsoil ! let soil liquid water equal to soil total water do i = 1, im if (islmsk(i) == 1) then slsoil(i,k) = smsoil(i,k) endif enddo enddo endif ! end if_lsm !!! !!! update surface diagnosis fields at the end of phys package !!! this change allows gocart to use filtered wind fields !!! if ( lgocart ) then call sfc_diag(im,pgr,gu0,gv0,gt0,gq0, & & tsea,qss,f10m,u10m,v10m,t2m,q2m,work3, & & evap,ffmm,ffhh,fm10,fh2) if (lssav) then do i = 1, im tmpmax(i) = max(tmpmax(i),t2m(i)) tmpmin(i) = min(tmpmin(i),t2m(i)) spfhmax(i) = max(spfhmax(i),q2m(i)) spfhmin(i) = min(spfhmin(i),q2m(i)) enddo endif endif ! --- ... total runoff is composed of drainage into water table and ! runoff at the surface and is accumulated in unit of meters if (lssav) then tem = dtf * 0.001 do i = 1, im runoff(i) = runoff(i) + (drain(i)+runof(i)) * tem srunoff(i) = srunoff(i) + runof(i) * tem enddo endif ! if (lprnt) write(0,*) ' after srunoff' ! --- ... xw: return updated ice thickness & concentration to global array do i = 1, im if (islmsk(i) == 2) then hice(i) = zice(i) fice(i) = cice(i) tisfc(i) = tice(i) else hice(i) = 0.0 fice(i) = 0.0 tisfc(i) = tsea(i) endif enddo ! --- ... return updated smsoil and stsoil to global arrays do k = 1, lsoil do i = 1, im smc(i,k) = smsoil(i,k) stc(i,k) = stsoil(i,k) slc(i,k) = slsoil(i,k) enddo enddo ! --- ... calculate column precipitable water "pwat" do i = 1, im pwat(i) = 0. ! rqtk(i) = 0. ! work2(i) = 1.0 / pgr(i) enddo do k = 1, levs do i = 1, im work1(i) = 0.0 enddo if (ncld > 0) then do ic = ntcw, ntcw+ncld-1 do i = 1, im work1(i) = work1(i) + gq0(i,k,ic) enddo enddo endif do i = 1, im ! if (lat == 45 .and. i == 1) write(1000+me,*)' pwatb=',pwat(1), ! &' kdt=',kdt,'del=',del(1,k),' gq0=',gq0(1,k,1),' work1=', ! &work1(1a,' k=',k pwat(i) = pwat(i) + del(i,k)*(gq0(i,k,1)+work1(i)) rqtk(i) = rqtk(i) + del(i,k)*(gq0(i,k,1)-qgrs(i,k,1)) enddo enddo do i = 1, im pwat(i) = pwat(i) * (1.0/con_g) ! rqtk(i) = rqtk(i) * work2(i) enddo ! ! if (lat == 45) write(1000+me,*)' pwat=',pwat(1),' kdt=',kdt ! if (lprnt) then ! write(0,*) ' endgt0=',gt0(ipr,:),' kdt=',kdt ! write(0,*) ' endgq0=',gq0(ipr,:,1),' kdt=',kdt ! endif deallocate (clw) if (do_shoc) then deallocate (qpl, qpi) endif if (.not. ras .or. .not. cscnv) then deallocate (cnvc, cnvw) endif ! deallocate (fscav, fswtr) ! ! if (lprnt) call mpi_quit(7) ! if (kdt >300 ) call mpi_quit(7) return !................................... end subroutine gbphys !-----------------------------------