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------------------------------- c include 'pars.h' c include 'parm.h' ! lq, lp, dx, dy, plon, plat c include 'coms.h' c include 'domain.h' ! slon, slat c include 'aghdr.h' ! IBYY, IBMM, IBDD, IBHH c include 'apraddcb.h' ! kountr integer lq, lp, dx, dy, plon, plat integer km integer slon, slat integer ibyy, ibmm, ibdd, ibhh integer khr real dt parameter( lq=1, lp=1, dx=0.5, dy=0.5 ) parameter( km=1 ) parameter( plon=lq, plat=lp ) parameter( ibyy=2001, ibmm=1, ibdd=1, ibhh=0 ) parameter( dt=450.0 ) parameter( khr=24, kfi=2 ) c real daysec ! seconds in one day = 86400.0 integer nht parameter(daysec=86400.0,nht=2280) ! c----------------------------------------------------------------------- integer nt ! count number of timesteps 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 ! day of SST file start (1st Wednesday) c character sst_source*3 ! REY or TMI c real dt_sst ! Fractal time interval in SST real aday0 ! real aday_model ! 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 SPdeS SST timing 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 c c vertical level from the top at 1 to the bottom at km c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c************************************************************* nhran=6 ! NCEP c nhran=12 !EC c sst_source = 'REY' ! for Reynolds, TMI for TMI daily c if ( sst_source .eq. 'REY' ) then nday_sst=7 ! weekly for fake Reynolds c else c nday_sst=1 ! daily c end if nhr_sst=nday_sst*24 dtan=real(nhran)*3600.0 c SPdeS SST timing c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c open(21,file='prep/sigma.d',form='formatted',status='unknown') c do k=1,km1 c read(21,5) kd, sig(k) c enddo c 5 format(2x,i2,4x,f10.6) c close(21) c do k=1,km c sig1(k)=0.5*(sig(k)+sig(k+1)) c enddo c c------------------------------------------------------------- c open(50,file='landsea/top.dat',form='unformatted', c + status='unknown') c read(50) os,lnd c close(50) c************************************************************* c Set model constants c------------------------------------------------------------- c c call cardin(ldbase,lsbase,lbdate,lbsec,kfi, c + finm,flname,restart,fbname) c call const(cb,bw,lq,lp,km,km1) c call param c------------------------------------------------------------- c Parameters used in Massflux convection scheme c------------------------------------------------------------- c call CUPARAM(km) c c-------------------------------- c to calculate Coriolis parameter c-------------------------------- c pi1=pi/180.0 c do 3 i=1,lp c y=slat+float(i-1)*dy c sni=sin(pi1*y) c f(i)=2.0*7.292e-5*sni c csf(i)=cos(pi1*y) c csfi(i)=1.0/csf(i) c tag(i)=sni*csfi(i)/6371.2e3 c 3 continue c++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c dti=1.0/dt inhr=int(3600.0/dt+0.5) c nttb=inhr*24*100 c 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 c if ( sst_source .eq. 'REY' ) then ibyy0sst=1990 ! fake Reynolds! ibmm0sst=1 ibdd0sst=3 ! Jan 3 is the first Wednesday in the data c else ! TMI c ibyy0sst=2001 c ibmm0sst=1 c 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 SPdeS SST timing c-------------------------------------------------------- c if(contn) then c open(40,file=restart,form='unformatted',status='old') c read(40) mtime,u,v,t,ps,qv,qc,qr,qi,qs,qg,omg,sst,ek,ep, c + zz0,el0,wst,sig1,sig,xlatent,sens,wsx,wsy,radt,glw, c + dir_vis,dif_vis,dir_nir,dif_nir,asdir,asdif,aldir, c + aldif,snow,t2m,rzsw,bht c close(40) c endif c nt=mtime*3600/int(dt+0.1) nhr=mtime nt0=nt nhr0=nhr c c SPdeS SST timing c lvlt=kday*(24/nhran)+mtime/nhran+1 c call initial(pso,uo,vo,to,qvo,os,lq,lp,km,dx,dy,81,lvlt, c + nhr,0,ibyy) c lvlt=1 c call initial(psn,un,vn,tn,qvn,os,lq,lp,km,dx,dy,81,1, c + nhr,1,ibyy) c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c if(.not.contn) then c do 43 i=1,lp c do 41 j=1,lq c ps(j,i)=pso(j,i) c asdir(j,i)=0.05 c asdif(j,i)=0.05 c aldir(j,i)=0.08 c aldif(j,i)=0.08 c asdirn(j,i)=0.05 c asdifn(j,i)=0.05 c aldirn(j,i)=0.08 c aldifn(j,i)=0.08 C for initialization of BATS related variables c xlatent(j,i) = 0.0 c sens(j,i) = 0.0 c xlatenta(j,i) = 0.0 c sensa(j,i) = 0.0 c wsx(j,i) = 0.1 c wsy(j,i) = 0.1 c wsxa(j,i) = 0.1 c wsya(j,i) = 0.1 c zz0(j,i)=1.0e-5 c el0(j,i)=1.0e6 c bht(j,i)=1000.0 c wst(j,i)=0.0 c snow(j,i)=1.e36 c t2m(j,i)=1.e36 c rzsw(j,i)=1.e36 c 41 continue c do 42 k=1,km c do 42 j=1,lq c u(j,i,k)=uo(j,i,k) c v(j,i,k)=vo(j,i,k) c t(j,i,k)=to(j,i,k) c qv(j,i,k)=qvo(j,i,k) c ek(j,i,k)=1.0e-4 c ep(j,i,k)=1.0e-6 c qc(j,i,k)=epsl c qr(j,i,k)=epsl c qi(j,i,k)=epsl c qs(j,i,k)=epsl c qg(j,i,k)=epsl c radt(j,i,k)=0.0 c radtn(j,i,k)=0.0 c 42 continue c 43 continue c endif c c Initialize radiation c ! c call inirad(lp, km, dy) c c++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c nyear90=ibyy-ibyy0sst kday_year=nyear90*365 kday_month=0 if(ibmm.gt.ibmm0sst) then do i=ibmm0sst,ibmm-1 kday_month=kday_month+mdat(i) enddo endif aday_model=1.0+3.5+real(ibhh+mtime)/24.0 ! HACK kday_model=int(aday_model) kdayt=kday_year+kday_month+kday_model+ibdd+nyear90/4 akss=real(kdayt)/real(nday_sst) ksst=int(akss) aday0=(akss-ksst)*nday_sst+(aday_model-kday_model) nday0=int(aday0) dt_sst=aday0/real(nday_sst) write(6,*) 'KSST, nday0,aday0,dt_sst ', + ksst,nday0,aday0,dt_sst c******************************************************** c c to read the SST fields (1990-p2000) c c-------------------------------------------------------- c if ( sst_source .eq. 'REY' ) then c call sst_set(ssto,lnd,lq,lp,dx,dy,79,ksst,0) c call sst_set(sstn,lnd,lq,lp,dx,dy,79,1 ,1) c else c call sst_set(ssto,lnd,lq,lp,dx,dy,39,ksst ,0) c kssti=ksst+1 c call sst_set(sstn,lnd,lq,lp,dx,dy,39,kssti,0) c end if c SPdeS in fact open a fake SST file that really contains dates c to test the temporal synchronization and interpolation. open(unit=39,file='fake_sst_date_file.dat',status='old', $ form='unformatted',access='direct',recl=4*1) c saved as 4 byte (32-bit) real read(unit=39,rec=ksst ) ssto kssti = ksst+1 read(unit=39,rec=kssti) sstn write(6,*) 'first two SST timestamps: ', ssto, sstn ! OK c if(.not.contn) then do 50 i=1,lp do 50 j=1,lq c if(lnd(j,i).eq.0) then sst(j,i)=ssto(j,i)+(sstn(j,i)-ssto(j,i))*dt_sst c else c sst(j,i)=t(j,i,km)*(1.0/sig1(km))**(rgas*gaa/grv) c endif 50 continue write(6,*) 'initial model interpolated time: ', sst ! OK c endif c c SPdeS SST timing c--------------------------------------------------- c initialize the model physics if not a restart run c--------------------------------------------------- c do 51 i=1,lp c do j=1,lq c rnl(j,i)=0.0 c rnc(j,i)=0.0 c olr(j,i)=0.0 c retswtop(j,i)=0.0 c clrolr(j,i)=0.0 c clrswtop(j,i)=0.0 c surlwd(j,i)=0.0 c clrsurlwd(j,i)=0.0 c gsw(j,i)=0.0 c clrswsur(j,i)=0.0 c topsw(j,i,1)=0.0 c topsw(j,i,2)=0.0 c topsw(j,i,3)=0.0 c topsw(j,i,4)=0.0 c psm(j,i)=0.0 c u10m(j,i)=0.0 c v10m(j,i)=0.0 c bhtm(j,i)=0.0 c sstm(j,i)=0.0 c cldtm(j,i)=0.0 c cldlm(j,i)=0.0 c sflxs(j,i)=0.0 c sflxl(j,i)=0.0 c sflxv(j,i)=0.0 c enddo c do k=1,km c do j=1,lq c um(j,i,k)=0.0 c vm(j,i,k)=0.0 c tm(j,i,k)=0.0 c qvm(j,i,k)=0.0 c qcm(j,i,k)=0.0 c qim(j,i,k)=0.0 c ekm(j,i,k)=0.0 c omgm(j,i,k)=0.0 c phim(j,i,k)=0.0 c dshtm(j,i,k)=0.0 c vdhtm(j,i,k)=0.0 c cvhtm(j,i,k)=0.0 c schtm(j,i,k)=0.0 c cdhtm(j,i,k)=0.0 c ttdm(j,i,k)=0.0 c rwradtm(j,i,k)=0.0 c swradtm(j,i,k)=0.0 c cldamtm(j,i,k)=0.0 c enddo c enddo c 51 continue c------------------------------------------------- c do j=1,lq c gglon(j)=slon+float(j-1)*dx c if(gglon(j).gt.180.0) gglon(j)=gglon(j)-360. c enddo c-------------------------------- c to manage the output filenames c-------------------------------- c nmn=100+mtime/khr c nmm=100+mtime/khr2 c------------------------------ c Initialization for BATS c------------------------------ c if(nt.eq.0) then c call zenith_drv(lq,lp,nt,0,gglon,azen,dfrac,dt) c call fakeini(nt,dt,lnd,contn,fbname,flname,azen,asdir, c + asdif,aldir,aldif,ldbase,lsbase,lbdate,lbsec,1) c endif c kfl=2 c if(contn) kfl=1 c call zenith_drv(lq,lp,nt,kountr,gglon,azenn,dfrac,dt) c call fakeini(nt,dt,lnd,contn,fbname,flname,azenn,asdirn, c + asdifn,aldirn,aldifn,ldbase,lsbase,lbdate,lbsec,kfl) c---------------------------------------- c Initialization for model physics c---------------------------------------- c if(.not.contn) then c qvb=qv c omg=0.0 c----------------------------------------------------------- c call phys(u,v,t,qv,qc,qr,qi,qs,qg,ek,ep,ps,rnl,rnc, c + zz0,el0,wst,sst,lnd,omg,qvb,lq,lp,km,km1,ut,vt,tt, c + radt,radtn,dsht,swradt,rwradt,cldamt,olr,retswtop, c + glw,glwn,xlatent,xlatenta,sens,sensa,wsx,wsy,nhr, c + dir_vis,dif_vis,dir_nir,dif_nir,dir_visn,dif_visn, c + dir_nirn,dif_nirn,asdir,asdif,aldir,aldif,asdirn, c + asdifn,aldirn,aldifn,snow,t2m,rzsw,gglon,azen,azenn, c + nmm,flname,nt0,clrolr,clrswtop,clrsurlwd,surlwd,gsw, c + clrswsur,vdht,cvht,cdht,scht,wsxa,wsya,u10,v10,bht, c + topsw) c write(6,*) 'Finishing initialization of model physics' c 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 c SPdeS SST timing c do 92 i=1,lp c do j=1,lq c pso(j,i)=psn(j,i) c enddo c do k=1,km c do j=1,lq c uo(j,i,k)=un(j,i,k) c vo(j,i,k)=vn(j,i,k) c to(j,i,k)=tn(j,i,k) c qvo(j,i,k)=qvn(j,i,k) c enddo c enddo c 92 continue c call initial(psn,un,vn,tn,qvn,os,lq,lp,km,dx,dy,81,1, c + nhr,1,ibyy) if(mod(nday0*24+nhr-nhr0,nhr_sst).eq.0) then do i=1,lp do j=1,lq ssto(j,i)=sstn(j,i) enddo 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 kssti=kssti+1 read(unit=39,rec=kssti) sstn endif c endif c c SPdeS SST timing c++++++++++++++++++++++++++++++++++++++++++++++++ c to write the model fields at every khr hours c------------------------------------------------ c if(mod(nhr,khr).eq.0) then c c nmn=nmn+1 c write(mmf,'(i3)') nmn c fname='/raid1/yqwang/result/'//finm//mmf//'.d' c fname='result/'//finm//mmf//'.d' c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c open(31,file=fname,form='unformatted',status='unknown') c write(31) nhr, c + um(stlon:enlon,stlat:enlat,:), c + vm(stlon:enlon,stlat:enlat,:), c + tm(stlon:enlon,stlat:enlat,:), c + psm(stlon:enlon,stlat:enlat), c + qvm(stlon:enlon,stlat:enlat,:), c + qcm(stlon:enlon,stlat:enlat,:), c + qim(stlon:enlon,stlat:enlat,:), c + omgm(stlon:enlon,stlat:enlat,:), c + phim(stlon:enlon,stlat:enlat,:), c + ekm(stlon:enlon,stlat:enlat,:), c + rnl(stlon:enlon,stlat:enlat), c + rnc(stlon:enlon,stlat:enlat), c + sstm(stlon:enlon,stlat:enlat), c + sig1, c + sig, c + sflxs(stlon:enlon,stlat:enlat), c + sflxl(stlon:enlon,stlat:enlat), c + sflxv(stlon:enlon,stlat:enlat), c + dshtm(stlon:enlon,stlat:enlat,:), c + swradtm(stlon:enlon,stlat:enlat,:), c + rwradtm(stlon:enlon,stlat:enlat,:), c + cldamtm(stlon:enlon,stlat:enlat,:), c + olr(stlon:enlon,stlat:enlat), c + retswtop(stlon:enlon,stlat:enlat), c + clrolr(stlon:enlon,stlat:enlat), c + clrswtop(stlon:enlon,stlat:enlat), c + clrsurlwd(stlon:enlon,stlat:enlat), c + surlwd(stlon:enlon,stlat:enlat), c + gsw(stlon:enlon,stlat:enlat), c + clrswsur(stlon:enlon,stlat:enlat), c + vdhtm(stlon:enlon,stlat:enlat,:), c + cvhtm(stlon:enlon,stlat:enlat,:), c + cdhtm(stlon:enlon,stlat:enlat,:), c + schtm(stlon:enlon,stlat:enlat,:), ! SPdeS changed to mean c + ttdm(stlon:enlon,stlat:enlat,:), c + u10m(stlon:enlon,stlat:enlat), c + v10m(stlon:enlon,stlat:enlat), c + bhtm(stlon:enlon,stlat:enlat), c + topsw(stlon:enlon,stlat:enlat,:), ! 4 bands? c + cldtm(stlon:enlon,stlat:enlat), c + cldlm(stlon:enlon,stlat:enlat) c close(31) c call flush(6) c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c to clear the khr-h accumulated rainfall for next khr-h use c------------------------------------------------------------ c do 54 i=1,lp c do j=1,lq c rnl(j,i)=0.0 c rnc(j,i)=0.0 c olr(j,i)=0.0 c retswtop(j,i)=0.0 c clrolr(j,i)=0.0 c clrswtop(j,i)=0.0 c surlwd(j,i)=0.0 c clrsurlwd(j,i)=0.0 c gsw(j,i)=0.0 c clrswsur(j,i)=0.0 c topsw(j,i,1)=0.0 c topsw(j,i,2)=0.0 c topsw(j,i,3)=0.0 c topsw(j,i,4)=0.0 c psm(j,i)=0.0 c u10m(j,i)=0.0 c v10m(j,i)=0.0 c bhtm(j,i)=0.0 c cldtm(j,i)=0.0 c cldlm(j,i)=0.0 c sstm(j,i)=0.0 c sflxs(j,i)=0.0 c sflxl(j,i)=0.0 c sflxv(j,i)=0.0 c enddo c do k=1,km c do j=1,lq c um(j,i,k)=0.0 c vm(j,i,k)=0.0 c tm(j,i,k)=0.0 c qvm(j,i,k)=0.0 c qcm(j,i,k)=0.0 c qim(j,i,k)=0.0 c ekm(j,i,k)=0.0 c omgm(j,i,k)=0.0 c phim(j,i,k)=0.0 c dshtm(j,i,k)=0.0 c vdhtm(j,i,k)=0.0 c cvhtm(j,i,k)=0.0 c schtm(j,i,k)=0.0 c cdhtm(j,i,k)=0.0 c ttdm(j,i,k)=0.0 c rwradtm(j,i,k)=0.0 c swradtm(j,i,k)=0.0 c cldamtm(j,i,k)=0.0 c enddo c enddo c 54 continue c c endif c c----------------------------------------------------------------- c To write the file for restart at time interval of khr2 in hours c----------------------------------------------------------------- c c if(mod(nhr,khr2).eq.0) then c c write(mmf,'(i3)') nmm c write(6,*) 'Writing restart file: nmm = ', nmm, ', mmf = ', mmf c rest_file='result/rest_'//finm//mmf c open(40,file=rest_file,form='unformatted',status='unknown') c write(40) nhr,u,v,t,ps,qv,qc,qr,qi,qs,qg,omg,sst,ek,ep, c + zz0,el0,wst,sig1,sig,xlatent,sens,wsx,wsy,radt,glw, c + dir_vis,dif_vis,dir_nir,dif_nir,asdir,asdif,aldir, c + aldif,snow,t2m,rzsw,bht c close(40) c c endif c c++++++++++++++++++++++++++++++++++++++++++++++++++++ c to quit the time integration c---------------------------------------------------- if(nhr.ge.kfi) goto 910 c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c 102 continue c c SPdeS SST timing c++++++++++++++++++++++++++++++++++++++++++++++++++ c integration core of the model c++++++++++++++++++++++++++++++++++++++++++++++++++ c do i=1,lp c do j=1,lq c c cldji=1.0 c do k=2,km c cldji=cldji*(1.0-cldamt(j,i,k)) c enddo c cldt(j,i)=1.0-cldji c c cldlji=1.0 c do k=17,km c cldlji=cldlji*(1.0-cldamt(j,i,k)) c enddo c cldl(j,i)=1.0-cldlji c enddo c enddo c++++++++++++++++++++++++++++++++++++++++++++++++++ c qvb=qv c ttd=t c-------------------------------------------------- c call inner(u,v,t,qv,qc,qr,qi,qs,qg,ek,ep,omg,ps,f, c + bw,csf,csfi,tag,lq,lp,km,km1,dx,dy,os,ut,vt,tt, c + uo,vo,to,qvo,pso,un,vn,tn,qvn,psn,phi,cb,dtan) c-------------------------------------------------------- c do i=1,lp c do k=1,km c do j=1,lq c ttd(j,i,k)=(t(j,i,k)-ttd(j,i,k))*dti-tt(j,i,k) c enddo c enddo c enddo c-------------------------------------------------------- c if(mod(nt,inhr).eq.0) then c nho=nt/inhr c c do k=1,km c asumu=0.0 c asumv=0.0 c asumt=0.0 c asumqv=0.0 c asumqc=0.0 c asumek=0.0 c asumsw=0.0 c asumlw=0.0 c asumch=0.0 c asumvh=0.0 c asumdh=0.0 c asumcv=0.0 c asumsc=0.0 c asumtt=0.0 c asumcl=0.0 cc do i=47,51 cc do j=117,121 c do i=81,85 c do j=161,165 c asumu=asumu+u(j,i,k) c asumv=asumv+v(j,i,k) c asumt=asumt+t(j,i,k) c asumqv=asumqv+qv(j,i,k) c asumqc=asumqc+qc(j,i,k) c asumek=asumek+ek(j,i,k) c asumsw=asumsw+swradt(j,i,k) c asumlw=asumlw+rwradt(j,i,k) c asumch=asumch+cdht(j,i,k) c asumdh=asumdh+dsht(j,i,k) c asumvh=asumvh+vdht(j,i,k) c asumcv=asumcv+cvht(j,i,k) c asumsc=asumsc+scht(j,i,k) c asumcl=asumcl+cldamt(j,i,k) c asumtt=asumtt+ttd(j,i,k) c enddo c enddo c uh(nho,k)=asumu/25.0 c vh(nho,k)=asumv/25.0 c th(nho,k)=asumt/25.0 c qvh(nho,k)=asumqv/25.0 c qch(nho,k)=asumqc/25.0 c ekh(nho,k)=asumek/25.0 c swh(nho,k)=asumsw/25.0 c lwh(nho,k)=asumlw/25.0 c cdhth(nho,k)=asumch/25.0 c vdhth(nho,k)=asumvh/25.0 c dshth(nho,k)=asumdh/25.0 c cvhth(nho,k)=asumcv/25.0 c schth(nho,k)=asumsc/25.0 c ttdh(nho,k)=asumtt/25.0 c cldh(nho,k)=asumcl/25.0 c enddo c asumps=0.0 c asumsh=0.0 c asumlh=0.0 c asumst=0.0 c asumbh=0.0 cc do i=47,51 cc do j=117,121 c do i=81,85 c do j=161,165 c asumps=asumps+ps(j,i) c asumsh=asumsh+sensa(j,i) c asumlh=asumlh+xlatenta(j,i) c asumst=asumst+sqrt(wsxa(j,i)*wsxa(j,i)+ c + wsya(j,i)*wsya(j,i)) c asumbh=asumbh+bht(j,i) c enddo c enddo c psh(nho)=asumps/25.0 c sflxsh(nho)=asumsh/25.0 c sflxlh(nho)=asumlh/25.0 c sflxvh(nho)=asumst/25.0 c bhth(nho)=asumbh/25.0 cc endif cc c do 55 i=1,lp c do j=1,lq c psm(j,i)=psm(j,i)+ps(j,i)*asum_sec c u10m(j,i)=u10m(j,i)+u10(j,i)*asum_sec c v10m(j,i)=v10m(j,i)+v10(j,i)*asum_sec c bhtm(j,i)=bhtm(j,i)+bht(j,i)*asum_sec c cldtm(j,i)=cldtm(j,i)+cldt(j,i)*asum_sec c cldlm(j,i)=cldlm(j,i)+cldl(j,i)*asum_sec c enddo c do k=1,km c do j=1,lq c um(j,i,k)=um(j,i,k)+u(j,i,k)*asum_sec c vm(j,i,k)=vm(j,i,k)+v(j,i,k)*asum_sec c tm(j,i,k)=tm(j,i,k)+t(j,i,k)*asum_sec c qvm(j,i,k)=qvm(j,i,k)+qv(j,i,k)*asum_sec c qcm(j,i,k)=qcm(j,i,k)+qc(j,i,k)*asum_sec c qim(j,i,k)=qim(j,i,k)+qi(j,i,k)*asum_sec c ekm(j,i,k)=ekm(j,i,k)+ek(j,i,k)*asum_sec c omgm(j,i,k)=omgm(j,i,k)+omg(j,i,k)*asum_sec c phim(j,i,k)=phim(j,i,k)+phi(j,i,k)*asum_sec c cdhtm(j,i,k)=cdhtm(j,i,k)+cdht(j,i,k)*asum_sec c vdhtm(j,i,k)=vdhtm(j,i,k)+vdht(j,i,k)*asum_sec c dshtm(j,i,k)=dshtm(j,i,k)+dsht(j,i,k)*asum_sec c cvhtm(j,i,k)=cvhtm(j,i,k)+cvht(j,i,k)*asum_sec c schtm(j,i,k)=schtm(j,i,k)+scht(j,i,k)*asum_sec c ttdm(j,i,k)=ttdm(j,i,k)+ttd(j,i,k)*asum_sec c swradtm(j,i,k)=swradtm(j,i,k)+swradt(j,i,k)*asum_sec c rwradtm(j,i,k)=rwradtm(j,i,k)+rwradt(j,i,k)*asum_sec c cldamtm(j,i,k)=cldamtm(j,i,k)+cldamt(j,i,k)*asum_sec c enddo c enddo c 55 continue cc--------------------------------------------------------- cc c do i=1,lp c do j=1,lq c if(lnd(j,i).eq.0) then c sflxs(j,i)=sflxs(j,i)+sensa(j,i)*asum_sec c sflxl(j,i)=sflxl(j,i)+xlatenta(j,i)*asum_sec c sflxv(j,i)=sflxv(j,i)+sqrt(wsxa(j,i)*wsxa(j,i)+wsya(j,i)* c + wsya(j,i))*asum_sec c else c sflxs(j,i)=sflxs(j,i)+sens(j,i)*asum_sec c sflxl(j,i)=sflxl(j,i)+xlatent(j,i)*asum_sec c sflxv(j,i)=sflxv(j,i)+sqrt(wsx(j,i)*wsx(j,i)+wsy(j,i)* c + wsy(j,i))*asum_sec c endif c enddo c enddo cc c----------------------------------------------------------- c call phys(u,v,t,qv,qc,qr,qi,qs,qg,ek,ep,ps,rnl,rnc, c + zz0,el0,wst,sst,lnd,omg,qvb,lq,lp,km,km1,ut,vt,tt, c + radt,radtn,dsht,swradt,rwradt,cldamt,olr,retswtop, c + glw,glwn,xlatent,xlatenta,sens,sensa,wsx,wsy,nhr, c + dir_vis,dif_vis,dir_nir,dif_nir,dir_visn,dif_visn, c + dir_nirn,dif_nirn,asdir,asdif,aldir,aldif,asdirn, c + asdifn,aldirn,aldifn,snow,t2m,rzsw,gglon,azen,azenn, c + nmm,flname,nt0,clrolr,clrswtop,clrsurlwd,surlwd,gsw, c + clrswsur,vdht,cvht,cdht,scht,wsxa,wsya,u10,v10,bht, c + topsw) c------------------------------------------------------------ dt_sst=amod(aday0*daysec+real(nt-nt0)*dt, !nday0 --> aday0 SPdeS + real(nday_sst)*daysec)/(real(nday_sst)*daysec) write(6,*) 'dt_sst = ', dt_sst if(dt_sst.lt.1.0e-20) dt_sst=1.0 do 203 i=1,lp do 203 j=1,lq c if(lnd(j,i).eq.0) then sst(j,i)=ssto(j,i)+(sstn(j,i)-ssto(j,i))*dt_sst c endif write(6,*) 'interp. time = ', sst sstm(j,i)=sstm(j,i)+sst(j,i)*asum_sec 203 continue 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 c SPdeS SST timing c diurn=finm//'_diurn.dat' c open(31,file=diurn,form='unformatted',status='unknown') c write(31) uh,vh,th,psh,qvh,qch,ekh,swh,lwh,cldh, c + cdhth,cvhth,schth,dshth,vdhth,ttdh,sflxsh,sflxlh, c + sflxvh,bhth,sig1,sig c close(31) c+++++++++++++++++++++++++++++++ c End the time integration c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c stop c end