C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c THIS MODEL BELONGS TO YUQING WANG, ANY ONE WHO WANTS TO USE THE MODEL% c SHOULD BE REQUEST TO HIM AND GET PERMISSION FROM HIM. IF ANY ERRORS % c ARE FOUND IN THE MODEL CODE PLEASE REPORT TO THE OWNER OF THE MODEL. % c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c*********************************************************************** c The model is designed for use in regional climate modeling using * c a multi-level hydrostatic primitive equation model with sigma as * c the vertical coordinate in sperical coordinates * c*********************************************************************** c c Yuqing Wang designed and coded the model (2001) c c Omer L. Sen helped in the implementation of the BATS land surface model c c Zhian Sun helped in the implementation of the radiation scheme c c======================================================================= c The main features of the model are as follows: c======================================================================= c c 1. Fourth-order conservative spacial finite-differencing scheme c 2. Combined Matsuno scheme and leapfrog scheme for time integration c 3. A modified Monin-Obkuhov similarity scheme for surface layer c 4. Either a level-1.5 or level-2 turbulence closure scheme for c vertical turbulent mixing above the surface layer. c 5. Biosphere-Atmosphere Transfer Scheme for land surface processes c 6. An advanced radiation scheme (Edward and Slingo, 1996) c c Details can be found from the BMRC Report by Yuqing Wang 1999. c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c program IPRC_RCMS c ================= C----------------------------------------------------------------------- implicit none C------------------------------Parameters------------------------------- include 'pars.h' include 'parm.h' ! lq, lp, dx, dy, plon, plat include 'coms.h' include 'domain.h' ! slon, slat include 'aghdr.h' ! IBYY, IBMM, IBDD, IBHH include 'apraddcb.h' ! kountr c real daysec ! seconds in one day = 86400.0 integer nht parameter(daysec=86400.0,nht=2280) ! c----------------------------------------------------------------------- integer ntime ! integer nday0 ! integer nmn ! Number for output writing integer nmm ! Number for restart output writing integer ksst ! timestamp for first sst read integer kssti ! incremented time stamp for sst read integer nyear90 ! integer kday_year ! integer kday_month ! integer kdayt ! integer kday_model ! integer nss ! integer nhr0 ! integer nt0 ! Restart time step integer lvlt ! integer j,i,k ! Loop parameters integer kd ! integer inhr ! integer mtime ! Time in hour integer nhr ! integer nho ! integer nhran ! Interval in hours for the intput driving fields integer nday_sst ! Interval in days for the driving SST integer nhr_sst ! Interval in hours for the driving SST integer kday ! integer kfi ! Total integration time in hour integer kfl ! integer ibyy0 ! year of reanalysis start integer ibmm0 ! month of reanalysis start integer ibyy0sst ! year of SST file start integer ibmm0sst ! month of SST file start integer ibdd0sst ! date of SST file start c character sst_source*3 ! REY or TMI c new work variables for SST timing SPdeS integer JD ! function JD for Julian date integer jdmodel0 ! Julian date of model start integration integer jdsst0 ! JD of SST data start real sstoffset ! diff between model start and SST start (const. days) real modeltime ! current time from model start (days) real ssttime ! current time from SST start (days) real sstotime ! time of older SST read real dayfracmodel0 ! fractional days of model start from yymmdd 00Z real dayfracsst0 ! frac. days of SST start from yymmdd 00Z real dt_sst ! Fractal time interval in SST real dt_ssto ! old dt_sst, for determining when to reread SST c real aday0 ! c real aday_model ! c real akss ! real dtan ! Interval in seconds for the input driving fields real dti ! 1.0/dt real dhx ! real dhy ! real pi1 ! real y ! real sni ! real asum_sec c real u(lq,lp,km) ! Zonal wind (m/s) real v(lq,lp,km) ! Meridional wind (m/s) real t(lq,lp,km) ! Air temperature (K) real ps(lq,lp) ! Surface Pressure (Pa) real f(lp) ! real omg(lq,lp,km) ! real qv(lq,lp,km) ! Water vapor mixing ratio (kg/kg) real zz0(lq,lp) ! Roughness length for momentum (m) real wst(lq,lp) ! real cb(lq,lp) ! real rnl(lq,lp) ! Large scale precipitation (mm/hr) real rnc(lq,lp) ! Convective precipitation (mm/hr) real el0(lq,lp) ! real qvb(lq,lp,km) ! previous time step water vapor mixing ratio real ek(lq,lp,km) ! real ep(lq,lp,km) ! real bw(lq,lp) ! real os(lq,lp) ! Topographic height (m) real sst(lq,lp) ! Sea surface temperature (K) real csf(lp) ! real csfi(lp) ! real tag(lp) ! real ut(lq,lp,km) ! real vt(lq,lp,km) ! real tt(lq,lp,km) ! real qc(lq,lp,km) ! real qr(lq,lp,km) ! real qi(lq,lp,km) ! real qs(lq,lp,km) ! real qg(lq,lp,km) ! real phi(lq,lp,km) ! geopotential height real gglon(lq) ! Longitudes c c Land surface model c real glw(lq,lp) ! Downward longwave radiation (W/m2) real glwn(lq,lp) ! Downward longwave radiation (W/m2) real dir_vis(lq,lp) ! Direct, visible radiation (W/m2) real dif_vis(lq,lp) ! Diffuse, visible radiation (W/m2) real dir_nir(lq,lp) ! Direct, NIR radiation (W/m2) real dif_nir(lq,lp) ! Diffuse, NIR radiation (W/m2) real dir_visn(lq,lp) ! Direct, visible radiation (W/m2) real dif_visn(lq,lp) ! Diffuse, visible radiation (W/m2) real dir_nirn(lq,lp) ! Direct, NIR radiation (W/m2) real dif_nirn(lq,lp) ! Diffuse, NIR radiation (W/m2) real xlatenta(lq,lp) ! Latent heat flux from "vdiff.f" (W/m2) real sensa(lq,lp) ! Sensible heat flux from "vdiff.f" (W/m2) real xlatent(lq,lp) ! Latent heat flux from BATS (W/m2) real sens(lq,lp) ! Sensible heat flux from BATS (W/m2) real wsx(lq,lp) ! Zonal wind stress real wsy(lq,lp) ! Meridional wind stress real wsxa(lq,lp) ! Zonal wind stress real wsya(lq,lp) ! Meridional wind stress real snow(lq,lp) ! Water equivalent snow (mm) real t2m(lq,lp) ! 2m air temperature (K) real rzsw(lq,lp) ! Root zone soil moisture real azen(lq,lp) ! Cosine of the solar zenith angle at time=t real dfrac(lq,lp) ! Day fraction real azenn(lq,lp) ! Cosine of the solar zenith angle at time=t+DT (radiation) integer lnd(lq,lp) ! Land/ocean index c integer ldbase ! Base day of run integer lsbase ! Base seconds of base day integer lbdate ! Base date of run (yymmdd format) integer lbsec ! Base seconds of base date c c Radiation and clouds c real asdir(lq,lp) ! Direct, visible albedo at time=t real asdif(lq,lp) ! Diffuse, visible albedo real aldir(lq,lp) ! Direct, NIR albedo real aldif(lq,lp) ! Diffuse, NIR albedo at time=t+DT (radiation) real asdirn(lq,lp) ! Direct, visible albedo real asdifn(lq,lp) ! Diffuse, visible albedo real aldirn(lq,lp) ! Direct, NIR albedo real aldifn(lq,lp) ! Diffuse, NIR albedo real dsht(lq,lp,km) ! dissipative heating rate real vdht(lq,lp,km) ! Vertical diffusion of T real cvht(lq,lp,km) ! Convective heating rate real scht(lq,lp,km) ! Shallow convective heating rate real cdht(lq,lp,km) ! Condensational heating rate real ttd(lq,lp,km) ! real rwradt(lq,lp,km) ! longwave heating rate real swradt(lq,lp,km) ! shortwave heating rate real retswtop(lq,lp) ! net shortwave flux at the model top real olr(lq,lp) ! outgoing longwave radiation real clrolr(lq,lp) ! Clear sky OLR real clrswtop(lq,lp) ! Clear sky net SW at the model top real surlwd(lq,lp) ! Downward LW radiation at the surface real clrsurlwd(lq,lp) ! Clear sky downward LW radiation at the surface real gsw(lq,lp) ! SW radiation at the surface real clrswsur(lq,lp) ! Clear sky SW radiation at the surface real topsw(lq,lp,4) ! SW fluxes at the model top real cldamt(lq,lp,km) ! cloud amount real cldt(lq,lp) ! Total cloudiness real cldtm(lq,lp) ! daily mean total cloudiness real cldl(lq,lp) ! low-level cloudiness real cldlm(lq,lp) ! low-level cloudiness real cldji, cldlji real radt(lq,lp,km) ! Net atmospheric radiation heating rate at time=t real radtn(lq,lp,km) ! Net atmospheric radiation heating rate at time=t+DT real um(lq,lp,km) ! real vm(lq,lp,km) ! real tm(lq,lp,km) ! real qvm(lq,lp,km) ! real qcm(lq,lp,km) ! real qim(lq,lp,km) ! real ekm(lq,lp,km) ! real omgm(lq,lp,km) ! real phim(lq,lp,km) ! real cdhtm(lq,lp,km) ! Condensational heating rate real cvhtm(lq,lp,km) ! Convective heating rate real schtm(lq,lp,km) ! Shallow convective heating rate real vdhtm(lq,lp,km) ! vertical mixing heating rate real dshtm(lq,lp,km) ! dissipative heating rate real ttdm(lq,lp,km) ! real rwradtm(lq,lp,km) ! longwave heating rate real swradtm(lq,lp,km) ! shortwave heating rate real cldamtm(lq,lp,km) ! cloud amount real psm(lq,lp) c c Lateral boundary update (two consequative to make interpolation) c real pso(lq,lp) ! Surface pressure (Pa) real uo(lq,lp,km) ! Zonal wind (m/s) real vo(lq,lp,km) ! Meridional wind (m/s) real to(lq,lp,km) ! Air temperature (K) real qvo(lq,lp,km) ! Water vapor mixing ratio (kg/kg) real psn(lq,lp) ! Surface pressure (Pa) real un(lq,lp,km) ! Zonal wind (m/s) real vn(lq,lp,km) ! Meridional wind (m/s) real tn(lq,lp,km) ! Air temperature (K) real qvn(lq,lp,km) ! Water vapor mixing ratio (kg/kg) real ssto(lq,lp) ! Sea surface temperature (K) real sstn(lq,lp) ! Sea surface temperature (K) c real coszro(lq) ! Cosine of the solar zenith angle real dyfrac(lq) ! Fraction of the day c real sflxs(lq,lp) ! Surface sensible heat flux real sflxl(lq,lp) ! Surface latent heat flux real sflxv(lq,lp) ! Surface wind stress real sstm(lq,lp) ! Surface temperature c real uh(nht,km) real vh(nht,km) real th(nht,km) real qvh(nht,km) real qch(nht,km) real ekh(nht,km) real swh(nht,km) real lwh(nht,km) real cldh(nht,km) real cdhth(nht,km) real vdhth(nht,km) real dshth(nht,km) real cvhth(nht,km) real schth(nht,km) real ttdh(nht,km) real psh(nht) real sflxsh(nht) real sflxlh(nht) real sflxvh(nht) real bhth(nht) ! Mixed-layer depth real u10(lq,lp) real v10(lq,lp) real bht(lq,lp) real u10m(lq,lp) real v10m(lq,lp) real bhtm(lq,lp) real asumu,asumv,asumt,asumqv,asumqc,asumek,asumlw,asumsw real asumch,asumtt,asumvh,asumps,asumsh,asumlh,asumst real asumcv,asumsc,asumdh,asumcl,asumbh c integer mdat(12) data mdat/31,28,31,30,31,30,31,31,30,31,30,31/ c character finm*5,fname*35,mmf*3,restart*40 character flname*5,fbname*15,rest_file*39,diurn*15 c c coordinates of the a subset of data to output in the .d files c (SPdeS 19 Nov 2001) c integer, parameter :: stlon=21,enlon=221 c integer, parameter :: stlat=11,enlat=131 integer, parameter :: stlon=1,enlon=lq integer, parameter :: stlat=1,enlat=lp c c vertical level from the top at 1 to the bottom at km c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c************************************************************* nhran=6 ! NCEP c nhran=12 !EC sst_source = 'TMI' ! REY for Reynolds, TMI for TMI daily if ( sst_source .eq. 'REY' ) then nday_sst=7 ! weekly else nday_sst=1 ! daily end if nhr_sst=nday_sst*24 dtan=real(nhran)*3600.0 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ open(21,file='prep/sigma.d',form='formatted',status='unknown') do k=1,km1 read(21,5) kd, sig(k) enddo 5 format(2x,i2,4x,f10.6) close(21) do k=1,km sig1(k)=0.5*(sig(k)+sig(k+1)) enddo c c------------------------------------------------------------- open(50,file='landsea/top.dat',form='unformatted', + status='unknown') read(50) os,lnd close(50) c************************************************************* c Set model constants c------------------------------------------------------------- c call cardin(ldbase,lsbase,lbdate,lbsec,kfi, + finm,flname,restart,fbname) call const(cb,bw,lq,lp,km,km1) call param c------------------------------------------------------------- c Parameters used in Massflux convection scheme c------------------------------------------------------------- call CUPARAM(km) c c-------------------------------- c to calculate Coriolis parameter c-------------------------------- pi1=pi/180.0 do 3 i=1,lp y=slat+float(i-1)*dy sni=sin(pi1*y) f(i)=2.0*7.292e-5*sni csf(i)=cos(pi1*y) csfi(i)=1.0/csf(i) tag(i)=sni*csfi(i)/6371.2e3 3 continue c++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c dti=1.0/dt inhr=int(3600.0/dt+0.5) c nttb=inhr*24*100 nttb=inhr c c-------------------------------------------------------- c to find the julian day to start this model simulation c-------------------------------------------------------- kday=0 ibmm0=1 if(ibyy.eq.1991) ibmm0=4 if(ibyy.eq.2001) then ibyy0 =2001 ibmm0 =1 endif dayfracmodel0=real(ibhh)/24.0 if ( sst_source .eq. 'REY' ) then ibyy0sst=1990 ibmm0sst=1 ibdd0sst=3 dayfracsst0=0.5 ! corresponds to 3 Jan 12:00 else ! TMI ibyy0sst=2001 ibmm0sst=1 ibdd0sst=1 dayfracsst0=0.5 ! SST data centered at noon of Jan 1 end if if(ibmm.gt.ibmm0) then do i=ibmm0,ibmm-1 kday=kday+mdat(i) enddo endif kday=kday+ibdd-1 mtime=0 c-------------------------------------------------------- if(contn) then open(40,file=restart,form='unformatted',status='old') read(40) mtime,u,v,t,ps,qv,qc,qr,qi,qs,qg,omg,sst,ek,ep, + zz0,el0,wst,sig1,sig,xlatent,sens,wsx,wsy,radt,glw, + dir_vis,dif_vis,dir_nir,dif_nir,asdir,asdif,aldir, + aldif,snow,t2m,rzsw,bht close(40) endif c nt=mtime*3600/int(dt+0.1) nhr=mtime nt0=nt nhr0=nhr write(6,*)'After restart: nt,nhr,mtime=',nt,nhr,mtime c lvlt=kday*(24/nhran)+mtime/nhran+1 call initial(pso,uo,vo,to,qvo,os,lq,lp,km,dx,dy,81,lvlt, + nhr,0,ibyy) lvlt=1 call initial(psn,un,vn,tn,qvn,os,lq,lp,km,dx,dy,81,1, + nhr,1,ibyy) c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ if(.not.contn) then do 43 i=1,lp do 41 j=1,lq ps(j,i)=pso(j,i) asdir(j,i)=0.05 asdif(j,i)=0.05 aldir(j,i)=0.08 aldif(j,i)=0.08 asdirn(j,i)=0.05 asdifn(j,i)=0.05 aldirn(j,i)=0.08 aldifn(j,i)=0.08 C for initialization of BATS related variables xlatent(j,i) = 0.0 sens(j,i) = 0.0 xlatenta(j,i) = 0.0 sensa(j,i) = 0.0 wsx(j,i) = 0.1 wsy(j,i) = 0.1 wsxa(j,i) = 0.1 wsya(j,i) = 0.1 zz0(j,i)=1.0e-5 el0(j,i)=1.0e6 bht(j,i)=1000.0 wst(j,i)=0.0 snow(j,i)=1.e36 t2m(j,i)=1.e36 rzsw(j,i)=1.e36 41 continue do 42 k=1,km do 42 j=1,lq u(j,i,k)=uo(j,i,k) v(j,i,k)=vo(j,i,k) t(j,i,k)=to(j,i,k) qv(j,i,k)=qvo(j,i,k) ek(j,i,k)=1.0e-4 ep(j,i,k)=1.0e-6 qc(j,i,k)=epsl qr(j,i,k)=epsl qi(j,i,k)=epsl qs(j,i,k)=epsl qg(j,i,k)=epsl radt(j,i,k)=0.0 radtn(j,i,k)=0.0 42 continue 43 continue endif c c Initialize radiation ! call inirad(lp, km, dy) c c++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c some constants for the sst timing jdmodel0= JD(ibyy ,ibmm ,ibdd ) jdsst0= JD(ibyy0sst,ibmm0sst,ibdd0sst) sstoffset= jdmodel0+dayfracmodel0 & -(jdsst0 +dayfracsst0 ) c check for model start before SST data if ( sstoffset .le. -nday_sst ) then write(6,*) 'Error: model start leads SST data by >= nday_sst.' stop elseif ( sstoffset .lt. 0.0 ) then write(6,*) 'Warning: model starts before SST data.' end if c the following times are modified each timestep nt modeltime= real(nt)*dt/daysec ssttime=modeltime+sstoffset c set ksst and dt_sst this way for the first time ksst= 1+max(0,int(ssttime/nday_sst)) sstotime= nday_sst*(ksst-1) dt_sst= (ssttime-sstotime)/nday_sst dt_ssto=dt_sst write(6,*) 'KSST,dt_sst ', + ksst,dt_sst ccc comment this out... c nyear90=ibyy-ibyy0sst c kday_year=nyear90*365 c kday_month=0 c if(ibmm.gt.ibmm0sst) then c do i=ibmm0sst,ibmm-1 c kday_month=kday_month+mdat(i) c enddo c endif c cc Reynolds SST elapsed time lags model time by 2.5 days, because this algorithm cc assumes the Reynolds SST data starts on 1 Jan 1990 midnight, but it cc actually starts on Jan 3, 1990 at noon (the first center of a time cc window). For TMI daily OI, the time interval is centered at noon, cc so day0sst=0.5 -SPdeS c c aday_model=-day0sst+real(ibhh+mtime)/24.0 c kday_model=int(aday_model) c kdayt=kday_year+kday_month+kday_model+ibdd+nyear90/4 c akss=real(kdayt)/real(nday_sst) + 1 ! SPdeS +1 changes to fort(1:n) subscript, rather than c(0:n-1) subscripts c ksst=int(akss) c aday0=(akss-ksst)*nday_sst+(aday_model-kday_model) c nday0=int(aday0) c dt_sst=aday0/real(nday_sst) c dt_ssto=dt_sst c write(6,*) 'KSST, nday0,aday0,dt_sst ', c + ksst,nday0,aday0,dt_sst ccc ...comment out to here c******************************************************** c c to read the SST fields (1990-p2000) c c-------------------------------------------------------- if ( sst_source .eq. 'REY' ) then call sst_set_rey(ssto,lnd,lq,lp,dx,dy,79,ksst,0) call sst_set_rey(sstn,lnd,lq,lp,dx,dy,79,1 ,1) else call sst_set_tmi(ssto,lnd,lq,lp,dx,dy,39,ksst ,0) kssti=ksst+1 call sst_set_tmi(sstn,lnd,lq,lp,dx,dy,39,kssti,0) end if c if(.not.contn) then do 50 i=1,lp do 50 j=1,lq if(lnd(j,i).eq.0) then sst(j,i)=ssto(j,i)+(sstn(j,i)-ssto(j,i))*dt_sst else sst(j,i)=t(j,i,km)*(1.0/sig1(km))**(rgas*gaa/grv) endif 50 continue endif c--------------------------------------------------- c initialize the model physics if not a restart run c--------------------------------------------------- do 51 i=1,lp do j=1,lq rnl(j,i)=0.0 rnc(j,i)=0.0 olr(j,i)=0.0 retswtop(j,i)=0.0 clrolr(j,i)=0.0 clrswtop(j,i)=0.0 surlwd(j,i)=0.0 clrsurlwd(j,i)=0.0 gsw(j,i)=0.0 clrswsur(j,i)=0.0 topsw(j,i,1)=0.0 topsw(j,i,2)=0.0 topsw(j,i,3)=0.0 topsw(j,i,4)=0.0 psm(j,i)=0.0 u10m(j,i)=0.0 v10m(j,i)=0.0 bhtm(j,i)=0.0 sstm(j,i)=0.0 cldtm(j,i)=0.0 cldlm(j,i)=0.0 sflxs(j,i)=0.0 sflxl(j,i)=0.0 sflxv(j,i)=0.0 enddo do k=1,km do j=1,lq um(j,i,k)=0.0 vm(j,i,k)=0.0 tm(j,i,k)=0.0 qvm(j,i,k)=0.0 qcm(j,i,k)=0.0 qim(j,i,k)=0.0 ekm(j,i,k)=0.0 omgm(j,i,k)=0.0 phim(j,i,k)=0.0 dshtm(j,i,k)=0.0 vdhtm(j,i,k)=0.0 cvhtm(j,i,k)=0.0 schtm(j,i,k)=0.0 cdhtm(j,i,k)=0.0 ttdm(j,i,k)=0.0 rwradtm(j,i,k)=0.0 swradtm(j,i,k)=0.0 cldamtm(j,i,k)=0.0 enddo enddo 51 continue c------------------------------------------------- do j=1,lq gglon(j)=slon+float(j-1)*dx if(gglon(j).gt.180.0) gglon(j)=gglon(j)-360. enddo c-------------------------------- c to manage the output filenames c-------------------------------- nmn=100+mtime/khr nmm=100+mtime/khr2 c------------------------------ c Initialization for BATS c------------------------------ if(nt.eq.0) then call zenith_drv(lq,lp,nt,0,gglon,azen,dfrac,dt) call fakeini(nt,dt,lnd,contn,fbname,flname,azen,asdir, + asdif,aldir,aldif,ldbase,lsbase,lbdate,lbsec,1) endif kfl=2 if(contn) kfl=1 call zenith_drv(lq,lp,nt,kountr,gglon,azenn,dfrac,dt) call fakeini(nt,dt,lnd,contn,fbname,flname,azenn,asdirn, + asdifn,aldirn,aldifn,ldbase,lsbase,lbdate,lbsec,kfl) c---------------------------------------- c Initialization for model physics c---------------------------------------- if(.not.contn) then qvb=qv omg=0.0 c----------------------------------------------------------- call phys(u,v,t,qv,qc,qr,qi,qs,qg,ek,ep,ps,rnl,rnc, + zz0,el0,wst,sst,lnd,omg,qvb,lq,lp,km,km1,ut,vt,tt, + radt,radtn,dsht,swradt,rwradt,cldamt,olr,retswtop, + glw,glwn,xlatent,xlatenta,sens,sensa,wsx,wsy,nhr, + dir_vis,dif_vis,dir_nir,dif_nir,dir_visn,dif_visn, + dir_nirn,dif_nirn,asdir,asdif,aldir,aldif,asdirn, + asdifn,aldirn,aldifn,snow,t2m,rzsw,gglon,azen,azenn, + nmm,flname,nt0,clrolr,clrswtop,clrsurlwd,surlwd,gsw, + clrswsur,vdht,cvht,cdht,scht,wsxa,wsya,u10,v10,bht, + topsw) write(6,*) 'Finishing initialization of model physics' endif c********************************************************************* c c Start time integration cycles c c==================================================================== c-------------------------------------------------------------------- c asum_sec=1.0/float(inhr*khr) write(6,*) 'asum_sec=',asum_sec c 100 nt=nt+1 c if((mod(nt-1,inhr).ne.0).or.(nt.eq.nt0+1)) goto 102 c nhr=(nt-1)/inhr c if(mod(nhr,nhran).eq.0.and.nhr.lt.kfi) then c do 92 i=1,lp do j=1,lq pso(j,i)=psn(j,i) enddo do k=1,km do j=1,lq uo(j,i,k)=un(j,i,k) vo(j,i,k)=vn(j,i,k) to(j,i,k)=tn(j,i,k) qvo(j,i,k)=qvn(j,i,k) enddo enddo 92 continue call initial(psn,un,vn,tn,qvn,os,lq,lp,km,dx,dy,81,1, + nhr,1,ibyy) cccccc Move this to below, s.t. a new SST record is read c whenever dt_sst wraps back to near zero. SST read should c not be conditional on the reanalysis read anyway. c SPdeS c c if(mod(nday0*24+nhr-nhr0,nhr_sst).eq.0) then c do i=1,lp c do j=1,lq c ssto(j,i)=sstn(j,i) c enddo c enddo c if ( sst_source .eq. 'REY' ) then c call sst_set(sstn,lnd,lq,lp,dx,dy,79,1 ,1) c else c kssti=kssti+1 c call sst_set(sstn,lnd,lq,lp,dx,dy,39,kssti,0) c end if c endif cccccc endif c c++++++++++++++++++++++++++++++++++++++++++++++++ c to write the model fields at every khr hours c------------------------------------------------ if(mod(nhr,khr).eq.0) then c nmn=nmn+1 write(mmf,'(i3)') nmn c fname='/raid1/yqwang/result/'//finm//mmf//'.d' fname='result/'//finm//mmf//'.d' c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ open(31,file=fname,form='unformatted',status='unknown') write(31) nhr, + um(stlon:enlon,stlat:enlat,:), + vm(stlon:enlon,stlat:enlat,:), + tm(stlon:enlon,stlat:enlat,:), + psm(stlon:enlon,stlat:enlat), + qvm(stlon:enlon,stlat:enlat,:), + qcm(stlon:enlon,stlat:enlat,:), + qim(stlon:enlon,stlat:enlat,:), + omgm(stlon:enlon,stlat:enlat,:), + phim(stlon:enlon,stlat:enlat,:), + ekm(stlon:enlon,stlat:enlat,:), + rnl(stlon:enlon,stlat:enlat), + rnc(stlon:enlon,stlat:enlat), + sstm(stlon:enlon,stlat:enlat), + sig1, + sig, + sflxs(stlon:enlon,stlat:enlat), + sflxl(stlon:enlon,stlat:enlat), + sflxv(stlon:enlon,stlat:enlat), + dshtm(stlon:enlon,stlat:enlat,:), + swradtm(stlon:enlon,stlat:enlat,:), + rwradtm(stlon:enlon,stlat:enlat,:), + cldamtm(stlon:enlon,stlat:enlat,:), + olr(stlon:enlon,stlat:enlat), + retswtop(stlon:enlon,stlat:enlat), + clrolr(stlon:enlon,stlat:enlat), + clrswtop(stlon:enlon,stlat:enlat), + clrsurlwd(stlon:enlon,stlat:enlat), + surlwd(stlon:enlon,stlat:enlat), + gsw(stlon:enlon,stlat:enlat), + clrswsur(stlon:enlon,stlat:enlat), + vdhtm(stlon:enlon,stlat:enlat,:), + cvhtm(stlon:enlon,stlat:enlat,:), + cdhtm(stlon:enlon,stlat:enlat,:), + schtm(stlon:enlon,stlat:enlat,:), ! SPdeS changed to mean + ttdm(stlon:enlon,stlat:enlat,:), + u10m(stlon:enlon,stlat:enlat), + v10m(stlon:enlon,stlat:enlat), + bhtm(stlon:enlon,stlat:enlat), + topsw(stlon:enlon,stlat:enlat,:), ! 4 bands? + cldtm(stlon:enlon,stlat:enlat), + cldlm(stlon:enlon,stlat:enlat) close(31) call flush(6) c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c to clear the khr-h accumulated rainfall for next khr-h use c------------------------------------------------------------ do 54 i=1,lp do j=1,lq rnl(j,i)=0.0 rnc(j,i)=0.0 olr(j,i)=0.0 retswtop(j,i)=0.0 clrolr(j,i)=0.0 clrswtop(j,i)=0.0 surlwd(j,i)=0.0 clrsurlwd(j,i)=0.0 gsw(j,i)=0.0 clrswsur(j,i)=0.0 topsw(j,i,1)=0.0 topsw(j,i,2)=0.0 topsw(j,i,3)=0.0 topsw(j,i,4)=0.0 psm(j,i)=0.0 u10m(j,i)=0.0 v10m(j,i)=0.0 bhtm(j,i)=0.0 cldtm(j,i)=0.0 cldlm(j,i)=0.0 sstm(j,i)=0.0 sflxs(j,i)=0.0 sflxl(j,i)=0.0 sflxv(j,i)=0.0 enddo do k=1,km do j=1,lq um(j,i,k)=0.0 vm(j,i,k)=0.0 tm(j,i,k)=0.0 qvm(j,i,k)=0.0 qcm(j,i,k)=0.0 qim(j,i,k)=0.0 ekm(j,i,k)=0.0 omgm(j,i,k)=0.0 phim(j,i,k)=0.0 dshtm(j,i,k)=0.0 vdhtm(j,i,k)=0.0 cvhtm(j,i,k)=0.0 schtm(j,i,k)=0.0 cdhtm(j,i,k)=0.0 ttdm(j,i,k)=0.0 rwradtm(j,i,k)=0.0 swradtm(j,i,k)=0.0 cldamtm(j,i,k)=0.0 enddo enddo 54 continue c endif c c----------------------------------------------------------------- c To write the file for restart at time interval of khr2 in hours c----------------------------------------------------------------- c if(mod(nhr,khr2).eq.0) then c write(mmf,'(i3)') nmm write(6,*) 'Writing restart file: nmm = ', nmm, ', mmf = ', mmf rest_file='result/rest_'//finm//mmf open(40,file=rest_file,form='unformatted',status='unknown') write(40) nhr,u,v,t,ps,qv,qc,qr,qi,qs,qg,omg,sst,ek,ep, + zz0,el0,wst,sig1,sig,xlatent,sens,wsx,wsy,radt,glw, + dir_vis,dif_vis,dir_nir,dif_nir,asdir,asdif,aldir, + aldif,snow,t2m,rzsw,bht close(40) c endif c c++++++++++++++++++++++++++++++++++++++++++++++++++++ c to quit the time integration c---------------------------------------------------- if(nhr.ge.kfi) goto 910 c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c 102 continue c SST read after NCEP read and possible output, before integration core modeltime= real(nt)*dt/daysec ssttime=modeltime+sstoffset c Cray compiler requires these type conversions: dt_sst=amod(ssttime,real(nday_sst))/real(nday_sst) c but the following is sufficient on Sun compiler c dt_sst=mod(ssttime,nday_sst)/nday_sst c Read a new SST record when dt_sst wraps around to near zero c update ssto, sstn if( dt_sst .lt. dt_ssto) then ssto=sstn if ( sst_source .eq. 'REY' ) then call sst_set_rey(sstn,lnd,lq,lp,dx,dy,79,1 ,1) else kssti=kssti+1 call sst_set_tmi(sstn,lnd,lq,lp,dx,dy,39,kssti,0) end if endif c Interpolate SST B.C. every timestep, update sst sst=ssto+(sstn-ssto)*dt_sst ! rely on compiler to vectorize write(6,*) sst, dt_sst ! test diagnostic sstm=sstm+sst*asum_sec dt_ssto=dt_sst c++++++++++++++++++++++++++++++++++++++++++++++++++ c integration core of the model c++++++++++++++++++++++++++++++++++++++++++++++++++ do i=1,lp do j=1,lq c cldji=1.0 do k=2,km cldji=cldji*(1.0-cldamt(j,i,k)) enddo cldt(j,i)=1.0-cldji c cldlji=1.0 do k=17,km cldlji=cldlji*(1.0-cldamt(j,i,k)) enddo cldl(j,i)=1.0-cldlji enddo enddo c++++++++++++++++++++++++++++++++++++++++++++++++++ qvb=qv ttd=t c-------------------------------------------------- call inner(u,v,t,qv,qc,qr,qi,qs,qg,ek,ep,omg,ps,f, + bw,csf,csfi,tag,lq,lp,km,km1,dx,dy,os,ut,vt,tt, + uo,vo,to,qvo,pso,un,vn,tn,qvn,psn,phi,cb,dtan) c-------------------------------------------------------- do i=1,lp do k=1,km do j=1,lq ttd(j,i,k)=(t(j,i,k)-ttd(j,i,k))*dti-tt(j,i,k) enddo enddo enddo c-------------------------------------------------------- c if(mod(nt,inhr).eq.0) then c nho=nt/inhr c do k=1,km asumu=0.0 asumv=0.0 asumt=0.0 asumqv=0.0 asumqc=0.0 asumek=0.0 asumsw=0.0 asumlw=0.0 asumch=0.0 asumvh=0.0 asumdh=0.0 asumcv=0.0 asumsc=0.0 asumtt=0.0 asumcl=0.0 c do i=47,51 c do j=117,121 do i=81,85 do j=161,165 asumu=asumu+u(j,i,k) asumv=asumv+v(j,i,k) asumt=asumt+t(j,i,k) asumqv=asumqv+qv(j,i,k) asumqc=asumqc+qc(j,i,k) asumek=asumek+ek(j,i,k) asumsw=asumsw+swradt(j,i,k) asumlw=asumlw+rwradt(j,i,k) asumch=asumch+cdht(j,i,k) asumdh=asumdh+dsht(j,i,k) asumvh=asumvh+vdht(j,i,k) asumcv=asumcv+cvht(j,i,k) asumsc=asumsc+scht(j,i,k) asumcl=asumcl+cldamt(j,i,k) asumtt=asumtt+ttd(j,i,k) enddo enddo uh(nho,k)=asumu/25.0 vh(nho,k)=asumv/25.0 th(nho,k)=asumt/25.0 qvh(nho,k)=asumqv/25.0 qch(nho,k)=asumqc/25.0 ekh(nho,k)=asumek/25.0 swh(nho,k)=asumsw/25.0 lwh(nho,k)=asumlw/25.0 cdhth(nho,k)=asumch/25.0 vdhth(nho,k)=asumvh/25.0 dshth(nho,k)=asumdh/25.0 cvhth(nho,k)=asumcv/25.0 schth(nho,k)=asumsc/25.0 ttdh(nho,k)=asumtt/25.0 cldh(nho,k)=asumcl/25.0 enddo asumps=0.0 asumsh=0.0 asumlh=0.0 asumst=0.0 asumbh=0.0 c do i=47,51 c do j=117,121 do i=81,85 do j=161,165 asumps=asumps+ps(j,i) asumsh=asumsh+sensa(j,i) asumlh=asumlh+xlatenta(j,i) asumst=asumst+sqrt(wsxa(j,i)*wsxa(j,i)+ + wsya(j,i)*wsya(j,i)) asumbh=asumbh+bht(j,i) enddo enddo psh(nho)=asumps/25.0 sflxsh(nho)=asumsh/25.0 sflxlh(nho)=asumlh/25.0 sflxvh(nho)=asumst/25.0 bhth(nho)=asumbh/25.0 c endif c do 55 i=1,lp do j=1,lq psm(j,i)=psm(j,i)+ps(j,i)*asum_sec u10m(j,i)=u10m(j,i)+u10(j,i)*asum_sec v10m(j,i)=v10m(j,i)+v10(j,i)*asum_sec bhtm(j,i)=bhtm(j,i)+bht(j,i)*asum_sec cldtm(j,i)=cldtm(j,i)+cldt(j,i)*asum_sec cldlm(j,i)=cldlm(j,i)+cldl(j,i)*asum_sec enddo do k=1,km do j=1,lq um(j,i,k)=um(j,i,k)+u(j,i,k)*asum_sec vm(j,i,k)=vm(j,i,k)+v(j,i,k)*asum_sec tm(j,i,k)=tm(j,i,k)+t(j,i,k)*asum_sec qvm(j,i,k)=qvm(j,i,k)+qv(j,i,k)*asum_sec qcm(j,i,k)=qcm(j,i,k)+qc(j,i,k)*asum_sec qim(j,i,k)=qim(j,i,k)+qi(j,i,k)*asum_sec ekm(j,i,k)=ekm(j,i,k)+ek(j,i,k)*asum_sec omgm(j,i,k)=omgm(j,i,k)+omg(j,i,k)*asum_sec phim(j,i,k)=phim(j,i,k)+phi(j,i,k)*asum_sec cdhtm(j,i,k)=cdhtm(j,i,k)+cdht(j,i,k)*asum_sec vdhtm(j,i,k)=vdhtm(j,i,k)+vdht(j,i,k)*asum_sec dshtm(j,i,k)=dshtm(j,i,k)+dsht(j,i,k)*asum_sec cvhtm(j,i,k)=cvhtm(j,i,k)+cvht(j,i,k)*asum_sec schtm(j,i,k)=schtm(j,i,k)+scht(j,i,k)*asum_sec ttdm(j,i,k)=ttdm(j,i,k)+ttd(j,i,k)*asum_sec swradtm(j,i,k)=swradtm(j,i,k)+swradt(j,i,k)*asum_sec rwradtm(j,i,k)=rwradtm(j,i,k)+rwradt(j,i,k)*asum_sec cldamtm(j,i,k)=cldamtm(j,i,k)+cldamt(j,i,k)*asum_sec enddo enddo 55 continue c--------------------------------------------------------- c do i=1,lp do j=1,lq if(lnd(j,i).eq.0) then sflxs(j,i)=sflxs(j,i)+sensa(j,i)*asum_sec sflxl(j,i)=sflxl(j,i)+xlatenta(j,i)*asum_sec sflxv(j,i)=sflxv(j,i)+sqrt(wsxa(j,i)*wsxa(j,i)+wsya(j,i)* + wsya(j,i))*asum_sec else sflxs(j,i)=sflxs(j,i)+sens(j,i)*asum_sec sflxl(j,i)=sflxl(j,i)+xlatent(j,i)*asum_sec sflxv(j,i)=sflxv(j,i)+sqrt(wsx(j,i)*wsx(j,i)+wsy(j,i)* + wsy(j,i))*asum_sec endif enddo enddo c c----------------------------------------------------------- call phys(u,v,t,qv,qc,qr,qi,qs,qg,ek,ep,ps,rnl,rnc, + zz0,el0,wst,sst,lnd,omg,qvb,lq,lp,km,km1,ut,vt,tt, + radt,radtn,dsht,swradt,rwradt,cldamt,olr,retswtop, + glw,glwn,xlatent,xlatenta,sens,sensa,wsx,wsy,nhr, + dir_vis,dif_vis,dir_nir,dif_nir,dir_visn,dif_visn, + dir_nirn,dif_nirn,asdir,asdif,aldir,aldif,asdirn, + asdifn,aldirn,aldifn,snow,t2m,rzsw,gglon,azen,azenn, + nmm,flname,nt0,clrolr,clrswtop,clrsurlwd,surlwd,gsw, + clrswsur,vdht,cvht,cdht,scht,wsxa,wsya,u10,v10,bht, + topsw) c------------------------------------------------------------ c--------------------------------------------- c--------------------------------------------- c End one full time step for all mesh domains c--------------------------------------------- c c------------------------------------------------- c------------------------------------------------- c Next time step c------------------------------------------------- c------------------------------------------------- goto 100 c c+++++++++++++++++++++++++ c 910 continue c diurn=finm//'_diurn.dat' open(31,file=diurn,form='unformatted',status='unknown') write(31) uh,vh,th,psh,qvh,qch,ekh,swh,lwh,cldh, + cdhth,cvhth,schth,dshth,vdhth,ttdh,sflxsh,sflxlh, + sflxvh,bhth,sig1,sig close(31) c+++++++++++++++++++++++++++++++ c End the time integration c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c stop c end