#include "rundeck_opts.h" CAOO Just to test CVS PROGRAM GISS_modelE,226 !@sum MAIN GISS modelE main time-stepping routine !@auth Original Development Team !@ver 1.0 (Based originally on B399) USE FILEMANAGER, only : openunit,closeunit USE TIMINGS, only : ntimemax,ntimeacc,timing,timestr USE PARAM USE MODEL_COM USE DOMAIN_DECOMP, ONLY : init_app,grid,AM_I_ROOT,pack_data USE DOMAIN_DECOMP, ONLY : ESMF_BCAST USE DYNAMICS USE RAD_COM, only : dimrad_sv USE RANDOM USE GETTIME_MOD #if (defined TRACERS_ON) || (defined TRACERS_OCEAN) USE TRACER_COM, only: mtrace #endif USE DIAG_COM, only : oa,monacc,koa,acc_period USE SOIL_DRV, only: daily_earth, ground_e USE SUBDAILY, only : nsubdd,init_subdd,get_subdd,reset_subdd USE DIAG_SERIAL, only : print_diags #ifdef USE_FVCORE USE FV_INTERFACE_MOD, only: fv_core USE FV_INTERFACE_MOD, only: Initialize USE FV_INTERFACE_MOD, only: Run USE FV_INTERFACE_MOD, only: Checkpoint USE FV_INTERFACE_MOD, only: Finalize USE FV_INTERFACE_MOD, only: Compute_Tendencies USE FV_INTERFACE_MOD, only: init_app_clock c$$$ USE MODEL_COM, only: clock USE ESMF_MOD, only: ESMF_Clock USE ESMF_CUSTOM_MOD, Only: vm => modelE_vm #endif USE ATMDYN, only : DYNAM,CALC_TROP,PGRAD_PBL,SDRAG & ,DISSIP,FILTER,CALC_AMPK, COMPUTE_DYNAM_AIJ_DIAGNOSTICS & ,COMPUTE_WSAVE, getTotalEnergy, addEnergyAsDiffuseHeat #ifdef TRACERS_ON & ,trdynam #endif USE ATMDYN_QDYNAM, only : QDYNAM use soil_drv, only : conserv_wtg, conserv_htg IMPLICIT NONE INTEGER K,M,MSTART,MNOW,MODD5D,months,ioerr,Ldate,istart INTEGER iu_VFLXO,iu_ACC,iu_RSF,iu_ODA INTEGER :: mpi_err, MDUM = 0 REAL*8, DIMENSION(NTIMEMAX) :: PERCENT REAL*8 DTIME,TOTALT , oa_glob(im,jm,koa) CHARACTER aDATE*14 CHARACTER*8 :: flg_go='___GO___' ! green light integer :: iflag=1 external sig_stop_model C**** Command line options LOGICAL :: qcrestart=.false. CHARACTER*32 :: ifile real :: lat_min=-90.,lat_max=90.,longt_min=0.,longt_max=360. integer :: tloopbegin, tloopend #ifdef USE_FVCORE Character(Len=*), Parameter :: fv_config = 'fv_config.rc' Type (FV_CORE) :: fv Type (ESMF_CLOCK) :: clock character(len=28) :: fv_fname, fv_dfname #endif integer :: L real*8 :: initialTotalEnergy, finalTotalEnergy ! tmp arrays real*8 w_ghy_j_2(jm),w_ghy_j_1(jm), w_lake_j_2(jm),w_lake_j_1(jm) real*8 h_ghy_j_2(jm),h_ghy_j_1(jm), h_lake_j_2(jm),h_lake_j_1(jm) call init_app(grid,im,jm,lm) call alloc_drv() C**** C**** Processing command line options C**** call read_options( qcrestart, ifile ) if ( qcrestart ) then call print_restart_info call stop_model("Terminated normally: printed restart info",-1) endif C**** C**** INITIALIZATIONS C**** CALL TIMER (MNOW,MDUM) #ifdef USE_FVCORE CALL INPUT (istart,ifile,clock) #else CALL INPUT (istart,ifile) #endif C**** Set run_status to "run in progress" if(istart > 0) call write_run_status("Run in progress...",1) C**** C**** If run is already done, just produce diagnostic printout C**** IF (Itime.GE.ItimeE.and.Kradia.le.0) then ! includes ISTART<1 case if (ItimeE.gt.0) then months=(Jyear-Jyear0)*JMperY + JMON-JMON0 call aPERIOD (JMON0,JYEAR0,months,1,0, acc_period,Ldate) end if call print_diags(0) if(istart < 1) then call stop_model ('Finished post-processing',-1) else call stop_model ('The run has already completed',13) end if ! no output files are affected END IF IF (AM_I_ROOT()) Then open(3,file='flagGoStop',form='FORMATTED',status='REPLACE') write (3,'(A8)') flg_go close (3) END IF call sys_signal( 15, sig_stop_model ) ! works only on single CPU MSTART=MNOW DO M=1,NTIMEACC MSTART= MSTART-TIMING(M) END DO C**** INITIALIZE TIME PARAMETERS NSTEP=(Itime-ItimeI)*NIdyn MODD5K=1000 CALL DAILY(.false.) ! not end_of_day CALL daily_RAD(.false.) if (istart.le.9) call reset_diag(0) if(Kradia==10) call daily_OCEAN(.false.) ! to test OCLIM if (Kradia.le.0) then CALL daily_EARTH(.false.) ! not end_of_day CALL daily_OCEAN(.false.) ! not end_of_day CALL CALC_AMPK(LS1-1) #if (defined TRACERS_ON) || (defined TRACERS_OCEAN) CALL daily_tracer(0) #endif if (kradia.le.0) CALL CHECKT ('INPUT ') end if CALL UPDTYPE C**** C**** Initialize FV dynamical core (ESMF component) if requested C**** #ifdef USE_FVCORE Call Initialize(fv, istart, vm, grid%esmf_grid, clock,fv_config) #endif if (AM_I_ROOT()) * WRITE (6,'(A,11X,A4,I5,A5,I3,A4,I3,6X,A,I4,I10)') * '0NASA/GISS Climate Model (re)started', * 'Year',JYEAR,aMON,JDATE,', Hr',JHOUR, * 'Internal clock: DTsrc-steps since 1/1/',Iyear1,ITIME CALL TIMER (MNOW,MELSE) C**** C**** Open and position output history files if needed C**** C**** Monthly files if (Kradia.ne.0 .and. Kradia<10) then write(aDATE(1:7),'(a3,I4.4)') aMON(1:3),Jyear if (Kradia.gt.0) aDATE(4:7)=' ' call openunit(trim('RAD'//aDATE(1:7)),iu_RAD,.true.,.false.) if (Kradia.lt.0) call io_POS(iu_RAD,Itime-1,2*dimrad_sv,Nrad) end if C**** Files for an accumulation period (1-12 months) write(aDATE(1:7),'(a3,I4.4)') aMON0(1:3),Jyear0 if (Kvflxo.ne.0) then call openunit('VFLXO'//aDATE(1:7),iu_VFLXO,.true.,.false.) call io_POS(iu_VFLXO,Itime,2*im*jm*koa,Nday) ! real*8-dim -> 2* end if C**** Initiallise file for sub-daily diagnostics, controlled by C**** space-separated string segments in SUBDD & SUBDD1 in the rundeck call init_subdd(aDATE) call sys_flush(6) C**** C**** MAIN LOOP C**** call gettime(tloopbegin) DO WHILE (Itime.lt.ItimeE) c$$$ call test_save(__LINE__, itime) C**** Every Ndisk Time Steps (DTsrc), starting with the first one, C**** write restart information alternately onto 2 disk files IF (MOD(Itime-ItimeI,Ndisk).eq.0) THEN CALL RFINAL (IRAND) call set_param( "IRAND", IRAND, 'o' ) iu_RSF=-1 IF (AM_I_ROOT()) * call openunit(rsf_file_name(KDISK),iu_RSF,.true.,.false.) call io_rsf(iu_RSF,Itime,iowrite,ioerr) IF (AM_I_ROOT()) call closeunit(iu_RSF) #ifdef USE_FVCORE fv_fname='fv.' ; write(fv_fname(4:4),'(i1)') kdisk fv_dfname='dfv.' ; write(fv_dfname(5:5),'(i1)') kdisk call Checkpoint(fv, clock, fv_fname, fv_dfname) #endif if (AM_I_ROOT()) * WRITE (6,'(A,I1,45X,A4,I5,A5,I3,A4,I3,A,I8)') * '0Restart file written on fort.',KDISK,'Year', * JYEAR,aMON,JDATE,', Hr',JHOUR,' Internal clock time:',ITIME KDISK=3-KDISK CALL TIMER (MNOW,MELSE) END IF C**** THINGS THAT GET DONE AT THE BEGINNING OF EVERY DAY IF (MOD(Itime,NDAY).eq.0) THEN C**** INITIALIZE SOME DIAG. ARRAYS AT THE BEGINNING OF SPECIFIED DAYS if (kradia.le.0) call daily_DIAG C**** THINGS THAT GET DONE AT THE BEGINNING OF EVERY MONTH IF ( JDAY.eq.1+JDendOfM(Jmon-1) ) then write(aDATE(1:7),'(a3,I4.4)') aMON(1:3),Jyear if (Kradia.ne.0 .and. Kradia<10) then if (Kradia.gt.0) aDATE(4:7)=' ' call closeunit( iu_RAD ) call openunit(trim('RAD'//aDATE(1:7)),iu_RAD,.true.,.false.) end if C**** THINGS THAT GET DONE AT THE BEGINNING OF EVERY ACC.PERIOD months=(Jyear-Jyear0)*JMperY + JMON-JMON0 if ( months.ge.NMONAV ) then call reset_DIAG(0) if (Kvflxo.ne.0) then call closeunit( iu_VFLXO ) call openunit('VFLXO'//aDATE(1:7),iu_VFLXO,.true.,.false.) end if C**** reset sub-daily diag files call reset_subdd(aDATE) end if ! beginning of acc.period END IF ! beginning of month END IF ! beginning of day C**** C**** INTEGRATE DYNAMIC TERMS (DIAGA AND DIAGB ARE CALLED FROM DYNAM) C**** if(Kradia>9) go to 100 ! to test daily/monthly procedures fast CALL CHECKT ('DYNAM0') if (kradia.le.0) then ! full model,kradia le 0 MODD5D=MOD(Itime-ItimeI,NDA5D) IF (MODD5D.EQ.0) IDACC(7)=IDACC(7)+1 IF (MODD5D.EQ.0) CALL DIAG5A (2,0) IF (MODD5D.EQ.0) CALL DIAGCA (1) PTOLD = P ! save for clouds C**** Initialize pressure for mass fluxes used by tracers and Q PS (:,:) = P(:,:) C**** Initialise total energy (J/m^2) initialTotalEnergy = getTotalEnergy() #ifndef USE_FVCORE CALL DYNAM() #else ! Using FV instead IF (MOD(Itime-ItimeI,NDAA).eq.0) CALL DIAGA0 call Run(fv, clock) CALL SDRAG (DTsrc) if (MOD(Itime-ItimeI,NDAA).eq.0) THEN call DIAGA call DIAGB call EPFLUX (U,V,T,P) endif #endif C**** This fix adjusts thermal energy to conserve total energy TE=KE+PE C**** Currently energy is put in uniformly weighted by mass finalTotalEnergy = getTotalEnergy() call addEnergyAsDiffuseHeat(finalTotalEnergy - initialTotalEnergy) call COMPUTE_DYNAM_AIJ_DIAGNOSTICS(PUA, PVA, DT) SD_CLOUDS(:,:,:) = CONV(:,:,:) call COMPUTE_WSAVE(wsave, sda, T, PK, PEDN, NIdyn) C**** Scale WM mixing ratios to conserve liquid water !$OMP PARALLEL DO PRIVATE (L) DO L=1,LS1-1 WM(:,:,L)=WM(:,:,L)* (PTOLD/P) END DO !$OMP END PARALLEL DO CALL QDYNAM ! Advection of Q by integrated fluxes CALL TIMER (MNOW,MDYN) #ifdef TRACERS_ON CALL TrDYNAM ! tracer dynamics CALL TIMER (MNOW,MTRACE) #endif C**** C**** Calculate tropopause level and pressure C**** CALL CALC_TROP C**** calculate some dynamic variables for the PBL CALL PGRAD_PBL C**** calculate zenith angle for current time step CALL CALC_ZENITH_ANGLE CALL CHECKT ('DYNAM ') CALL TIMER (MNOW,MSURF) IF (MODD5D.EQ.0) CALL DIAG5A (7,NIdyn) IF (MODD5D.EQ.0) CALL DIAGCA (2) IF (MOD(Itime,NDAY/2).eq.0) CALL DIAG7A C**** C**** INTEGRATE SOURCE TERMS C**** IDACC(1)=IDACC(1)+1 MODD5S=MOD(Itime-ItimeI,NDA5S) IF (MODD5S.EQ.0) IDACC(8)=IDACC(8)+1 IF (MODD5S.EQ.0.AND.MODD5D.NE.0) CALL DIAG5A (1,0) IF (MODD5S.EQ.0.AND.MODD5D.NE.0) CALL DIAGCA (1) C**** FIRST CALL MELT_SI SO THAT TOO SMALL ICE FRACTIONS ARE REMOVED C**** AND ICE FRACTION CAN THEN STAY CONSTANT UNTIL END OF TIMESTEP CALL MELT_SI CALL UPDTYPE CALL TIMER (MNOW,MSURF) C**** CONDENSATION, SUPER SATURATION AND MOIST CONVECTION CALL CONDSE CALL CHECKT ('CONDSE') CALL TIMER (MNOW,MCNDS) IF (MODD5S.EQ.0) CALL DIAG5A (9,NIdyn) IF (MODD5S.EQ.0) CALL DIAGCA (3) end if ! full model,kradia le 0 C**** RADIATION, SOLAR AND THERMAL MODRD=MOD(Itime-ItimeI,NRAD) if (kradia.le.0. or. MODRD.eq.0) then CALL RADIA if (kradia.le.0) CALL CHECKT ('RADIA ') end if CALL TIMER (MNOW,MRAD) if (kradia.le.0) then ! full model,kradia le 0 IF (MODD5S.EQ.0) CALL DIAG5A (11,NIdyn) C**** C**** SURFACE INTERACTION AND GROUND CALCULATION C**** C**** NOTE THAT FLUXES ARE APPLIED IN TOP-DOWN ORDER SO THAT THE C**** FLUXES FROM ONE MODULE CAN BE SUBSEQUENTLY APPLIED TO THAT BELOW C**** IF (MODD5S.EQ.0) CALL DIAGCA (4) C**** APPLY PRECIPITATION TO SEA/LAKE/LAND ICE CALL PRECIP_SI CALL PRECIP_LI C**** APPLY PRECIPITATION AND RUNOFF TO LAKES/OCEANS CALL PRECIP_LK CALL PRECIP_OC CALL TIMER (MNOW,MSURF) CALL CHECKT ('PRECIP') #ifdef TRACERS_ON C**** Calculate non-interactive tracer surface sources and sinks call set_tracer_2Dsource CALL TIMER (MNOW,MTRACE) #endif C**** CALCULATE SURFACE FLUXES AND EARTH CALL SURFCE CALL CHECKT ('SURFCE') CALL TIMER (MNOW,MSURF) IF (MODD5S.EQ.0) CALL DIAGCA (5) C**** CALCULATE ICE DYNAMICS CALL DYNSI C**** CALCULATE BASE ICE-OCEAN/LAKE FLUXES CALL UNDERICE C**** APPLY SURFACE/BASE FLUXES TO SEA/LAKE ICE CALL GROUND_SI C**** APPLY SURFACE FLUXES TO LAND ICE CALL GROUND_LI CALL CHECKT ('GRNDSI') C**** APPLY FLUXES TO LAKES AND DETERMINE ICE FORMATION CALL GROUND_LK CALL CHECKT ('GRNDLK') CALL TIMER (MNOW,MSURF) IF (MODD5S.EQ.0) CALL DIAGCA (6) C**** CALCULATE RIVER RUNOFF FROM LAKE MASS CALL RIVERF CALL GROUND_E ! diagnostic only - should be merged with EARTH C**** APPLY FLUXES TO OCEAN, DO OCEAN DYNAMICS AND CALC. ICE FORMATION CALL OCEANS CALL CHECKT ('OCEANS') C**** APPLY ICE FORMED IN THE OCEAN/LAKES TO ICE VARIABLES CALL FORM_SI CALL CHECKT ('FORMSI') C**** IF ATURB is used in rundeck then this is a dummy call C**** CALCULATE DRY CONVECTION ABOVE PBL CALL ATM_DIFFUS (2,LM-1,dtsrc) CALL CHECKT ('DRYCNV') CALL TIMER (MNOW,MSURF) IF (MODD5S.EQ.0) CALL DIAGCA (9) C**** ADVECT ICE CALL ADVSI CALL ADVSI_DIAG ! needed to update qflux model, dummy otherwise CALL CHECKT ('ADVSI ') C**** UPDATE DIAGNOSTIC TYPES CALL UPDTYPE C**** ADD DISSIPATED KE FROM SURFACE CALCULATION BACK AS LOCAL HEAT CALL DISSIP CALL CHECKT ('DISSIP') CALL TIMER (MNOW,MSURF) IF (MODD5S.EQ.0) CALL DIAGCA (7) IF (MODD5S.EQ.0) CALL DIAG5A (12,NIdyn) C**** SEA LEVEL PRESSURE FILTER IF (MFILTR.GT.0.AND.MOD(Itime-ItimeI,NFILTR).EQ.0) THEN IDACC(10)=IDACC(10)+1 IF (MODD5S.NE.0) CALL DIAG5A (1,0) CALL DIAGCA (1) CALL FILTER CALL CHECKT ('FILTER') CALL TIMER (MNOW,MDYN) CALL DIAG5A (14,NFILTR*NIdyn) CALL DIAGCA (8) END IF #ifdef TRACERS_ON C**** 3D Tracer sources and sinks C**** Tracer gravitational settling for aerosols CALL TRGRAV C**** Tracer radioactive decay (and possible source) CALL TDECAY C**** Calculate 3D tracers sources and sinks call tracer_3Dsource C**** Accumulate tracer distribution diagnostics CALL TRACEA CALL TIMER (MNOW,MTRACE) CALL CHECKT ('T3DSRC') #endif end if ! full model,kradia le 0 C**** C**** WRITE SUB-DAILY DIAGNOSTICS EVERY NSUBDD hours C**** if (Nsubdd.ne.0) then if (mod(Itime+1,Nsubdd).eq.0) call get_subdd end if #ifdef TRACERS_DUST call ahourly #endif C**** C**** UPDATE Internal MODEL TIME AND CALL DAILY IF REQUIRED C**** 100 Itime=Itime+1 ! DTsrc-steps since 1/1/Iyear1 Jhour=MOD(Itime*24/NDAY,24) ! Hour (0-23) Nstep=Nstep+NIdyn ! counts DT(dyn)-steps IF (MOD(Itime,NDAY).eq.0) THEN ! NEW DAY if (kradia.gt.0) then ! radiative forcing run CALL DAILY(.false.) if(Kradia<10) CALL daily_RAD(.true.) if(Kradia==10) CALL daily_OCEAN(.true.) ! to test OCLIM months=(Jyear-Jyear0)*JMperY + JMON-JMON0 else ! full model, kradia le 0 CALL DIAG5A (1,0) CALL DIAGCA (1) CALL DAILY(.true.) ! end_of_day CALL daily_RAD(.true.) months=(Jyear-Jyear0)*JMperY + JMON-JMON0 CALL TIMER (MNOW,MELSE) cddd call conserv_wtg(w_ghy_j_1) cddd call conserv_LKM(w_lake_j_1) cddd call conserv_htg(h_ghy_j_1) cddd call conserv_LKE(h_lake_j_1) call daily_LAKE call daily_EARTH(.true.) ! end_of_day cddd call conserv_wtg(w_ghy_j_2) cddd call conserv_LKM(w_lake_j_2) cddd call conserv_htg(h_ghy_j_2) cddd call conserv_LKE(h_lake_j_2) cddd print *,"GHY_LAKE_CONSERV:", cddd & sum( w_ghy_j_1+w_lake_j_1 ), cddd & sum( w_ghy_j_2-w_ghy_j_1+w_lake_j_2-w_lake_j_1 ), cddd & sum( h_ghy_j_1+h_lake_j_1 ), cddd & sum( h_ghy_j_2-h_ghy_j_1+h_lake_j_2-h_lake_j_1 ) cddd print *,"GHY_LAKE tot hgy, tot lake, d ghy, d lake", cddd & sum( w_ghy_j_1 ), sum( w_lake_j_1 ), cddd & sum( w_ghy_j_2-w_ghy_j_1 ), cddd & sum( w_lake_j_2-w_lake_j_1 ) cddd ! print *,"GHY_LAKE errors:", cddd ! & w_ghy_j_2-w_ghy_j_1+w_lake_j_2-w_lake_j_1 cddd print *,"GHY_LAKE_E tot hgy, tot lake, d ghy, d lake", cddd & sum( h_ghy_j_1 ), sum( h_lake_j_1 ), cddd & sum( h_ghy_j_2-h_ghy_j_1 ), cddd & sum( h_lake_j_2-h_lake_j_1 ) cddd print *,"GHY_LAKE errors:", cddd & h_ghy_j_2-h_ghy_j_1+h_lake_j_2-h_lake_j_1 call daily_OCEAN(.true.) ! end_of_day call daily_ICE call daily_LI #if (defined TRACERS_ON) || (defined TRACERS_OCEAN) call daily_tracer(1) CALL TIMER (MNOW,MTRACE) #endif CALL CHECKT ('DAILY ') CALL TIMER (MNOW,MSURF) CALL DIAG5A (16,NDAY*NIdyn) CALL DIAGCA (10) call sys_flush(6) end if ! kradia: full model (or rad.forcing run) CALL UPDTYPE END IF ! NEW DAY #ifdef USE_FVCORE Call Compute_Tendencies(fv) #endif if (kradia.le.0) then ! full model C**** C**** WRITE INFORMATION FOR OHT CALCULATION EVERY 24 HOURS C**** IF (Kvflxo.EQ.0.) OA(:,:,4:KOA)=0. ! to prepare for future saves IF (Kvflxo.NE.0.) THEN IF (MOD(Itime,NDAY).eq.0) THEN call pack_data (grid, OA, OA_glob) if (am_I_root()) * call WRITEI8 (iu_vflxo,Itime,OA_glob,im*jm*koa) C**** ZERO OUT INTEGRATED QUANTITIES OA(:,:,4:KOA)=0. ELSEIF (MOD(Itime,NDAY/2).eq.0) THEN call vflx_OCEAN END IF CALL TIMER (MNOW,MELSE) END IF C**** C**** CALL DIAGNOSTIC ROUTINES C**** IF (MOD(Itime-ItimeI,NDA4).EQ.0) CALL DIAG4A ! at hr 23 E-history C**** PRINT CURRENT DIAGNOSTICS (INCLUDING THE INITIAL CONDITIONS) IF (NIPRNT.GT.0) THEN acc_period='PARTIAL ' call print_diags(1) NIPRNT=NIPRNT-1 call set_param( "NIPRNT", NIPRNT, 'o' ) END IF end if ! full model ; kradia le 0 C**** THINGS TO DO BEFORE ZEROING OUT THE ACCUMULATING ARRAYS C**** done at the end of (selected) months IF (months.ge.NMONAV .and. ! next 2 conditions are rarely needed * JDAY.eq.1+JDendOfM(JMON-1) .and. MOD(Itime,NDAY).eq.0) THEN C**** PRINT DIAGNOSTIC TIME AVERAGED QUANTITIES call aPERIOD (JMON0,JYEAR0,months,1,0, aDATE(1:12),Ldate) acc_period=aDATE(1:12) WRITE (aDATE(8:14),'(A3,I4.4)') aMON(1:3),JYEAR if (kradia.le.0) call print_diags(0) C**** SAVE ONE OR BOTH PARTS OF THE FINAL RESTART DATA SET IF (KCOPY.GT.0) THEN C**** KCOPY > 0 : SAVE THE DIAGNOSTIC ACCUM ARRAYS IN SINGLE PRECISION monacc = 0 do k=JMON0,JMON0+NMONAV-1 m = k if(m.gt.12) m = m-12 monacc(m) = 1 end do If (AM_I_ROOT()) * call openunit(aDATE(1:7)//'.acc'//XLABEL(1:LRUNID),iu_ACC, * .true.,.false.) call io_rsf (iu_ACC,Itime,iowrite_single,ioerr) If (AM_I_ROOT()) call closeunit(iu_ACC) C**** KCOPY > 1 : ALSO SAVE THE RESTART INFORMATION IF (KCOPY.GT.1) THEN CALL RFINAL (IRAND) call set_param( "IRAND", IRAND, 'o' ) IF (AM_I_ROOT()) * call openunit('1'//aDATE(8:14)//'.rsf'//XLABEL(1:LRUNID) * ,iu_RSF,.true.,.false.) call io_rsf(iu_RSF,Itime,iowrite_mon,ioerr) IF (AM_I_ROOT()) call closeunit(iu_RSF) #ifdef USE_FVCORE fv_fname = '1'//aDATE(8:14)//'.fv'//XLABEL(1:LRUNID) fv_dfname = '1'//aDATE(8:14)//'.dfv'//XLABEL(1:LRUNID) call Checkpoint(fv, clock, fv_fname, fv_dfname) #endif END IF C**** KCOPY > 2 : ALSO SAVE THE OCEAN DATA TO INITIALIZE DEEP OCEAN RUNS IF (KCOPY.GT.2) THEN If (AM_I_ROOT()) * call openunit(aDATE(1:7)//'.oda'//XLABEL(1:LRUNID) * ,iu_ODA,.true.,.false.) call io_oda(iu_ODA,Itime,iowrite,ioerr) IF (AM_I_ROOT()) call closeunit(iu_ODA) END IF END IF C**** PRINT AND ZERO OUT THE TIMING NUMBERS CALL TIMER (MNOW,MDIAG) TOTALT=SUM(TIMING(1:NTIMEACC)) DO M=1,NTIMEACC PERCENT(M) = 100d0*TIMING(M)/(TOTALT+.00001) END DO c TOTALT=(MNOW-MSTART)/(60.*100.) ! wrong when clock rolls over TOTALT=TOTALT/(60.*100.) ! in minutes DTIME = NDAY*TOTALT/(Itime-Itime0) ! minutes/day WRITE (6,'(/A,F7.2,A,/(8(A13,F5.1/))//)') * '0TIME',DTIME,'(MINUTES) ',(TIMESTR(M),PERCENT(M),M=1,NTIMEACC) TIMING = 0 MSTART= MNOW END IF C**** CPU TIME FOR CALLING DIAGNOSTICS CALL TIMER (MNOW,MDIAG) C**** TEST FOR TERMINATION OF RUN ccc IF (MOD(Itime,Nssw).eq.0) then IF (AM_I_ROOT()) then flg_go = '__STOP__' ! stop if flagGoStop if missing iflag=0 open(3,file='flagGoStop',form='FORMATTED',status='OLD',err=210) read (3,'(A8)',end=210) flg_go close (3) 210 continue IF (flg_go .eq. '___GO___') iflag=1 call ESMF_BCAST( grid, iflag) else call ESMF_BCAST( grid, iflag) if (iflag .eq. 1) flg_go = '___GO___' if (iflag .eq. 0) flg_go = '__STOP__' end if endif IF (flg_go.ne.'___GO___' .or. stop_on) THEN C**** Flag to continue run has been turned off WRITE (6,'("0Flag to continue run has been turned off.")') EXIT END IF c$$$ call test_save(__LINE__, itime-1) END DO call gettime(tloopend) if (AM_I_ROOT()) * write(*,*) "Time spent in the main loop in seconds:", * .01*(tloopend-tloopbegin) C**** C**** END OF MAIN LOOP C**** C**** ALWAYS PRINT OUT RSF FILE WHEN EXITING CALL RFINAL (IRAND) call set_param( "IRAND", IRAND, 'o' ) IF (AM_I_ROOT()) * call openunit(rsf_file_name(KDISK),iu_RSF,.true.,.false.) call io_rsf(iu_RSF,Itime,iowrite,ioerr) IF (AM_I_ROOT()) call closeunit(iu_RSF) #ifdef USE_FVCORE fv_fname='fv.' ; write(fv_fname(4:4),'(i1)') kdisk fv_dfname='dfv.' ; write(fv_dfname(5:5),'(i1)') kdisk call Finalize(fv, clock, fv_fname, fv_dfname) #endif if (AM_I_ROOT()) then WRITE (6,'(A,I1,45X,A4,I5,A5,I3,A4,I3,A,I8)') * '0Restart file written on fort.',KDISK,'Year',JYEAR, * aMON,JDATE,', Hr',JHOUR,' Internal clock time:',ITIME end if C**** RUN TERMINATED BECAUSE IT REACHED TAUE (OR SS6 WAS TURNED ON) IF (AM_I_ROOT()) * WRITE (6,'(/////4(1X,33("****")/)//,A,I8 * ///4(1X,33("****")/))') * ' PROGRAM TERMINATED NORMALLY - Internal clock time:',ITIME IF (Itime.ge.ItimeE) CALL stop_model ( & 'Terminated normally (reached maximum time)',13) CALL stop_model ('Run stopped with sswE',12) ! voluntary stop END subroutine sig_stop_model,1 USE MODEL_COM, only : stop_on implicit none stop_on = .true. end subroutine sig_stop_model subroutine init_Model 1,36 !@sum This program reads most of parameters from the database (DB) !@+ get_param( "A", X ) reads parameter A into variable X !@+ if "A" is not in the database, it will generate an error !@+ message and stop !@+ sync_param( "B", Y ) reads parameter B into variable Y !@+ if "B" is not in the database, then Y is unchanged and its !@+ value is saved in the database as "B" (here sync = synchronize) USE MODEL_COM, only : JM,LM,NIPRNT,MFILTR,NFILTR,NRAD * ,NDASF,NDA4,NDA5S,NDA5K,NDA5D,NDAA,Kvflxo,kradia * ,NMONAV,Ndisk,Nssw,KCOPY,KOCEAN,NIsurf,iyear1 $ ,LS1,IRAND,ItimeI,PSTRAT,UOdrag $ ,X_SDRAG,C_SDRAG,LSDRAG,P_SDRAG,LPSDRAG,PP_SDRAG,ang_sdrag $ ,P_CSDRAG,CSDRAGL,Wc_Jdrag,COUPLED_CHEM,dt * ,DT_XUfilter,DT_XVfilter,DT_YVfilter,DT_YUfilter,QUVfilter & ,do_polefix,pednl00,pmidl00 USE DOMAIN_DECOMP, only: AM_I_ROOT USE PARAM implicit none INTEGER L,LCSDRAG C**** Rundeck parameters: call sync_param( "NMONAV", NMONAV ) call sync_param( "NIPRNT", NIPRNT ) call sync_param( "DT_XVfilter", DT_XVfilter ) call sync_param( "DT_XUfilter", DT_XUfilter ) call sync_param( "DT_YVfilter", DT_YVfilter ) call sync_param( "DT_YUfilter", DT_YUfilter ) call sync_param( "MFILTR", MFILTR ) call sync_param( "X_SDRAG", X_SDRAG, 2 ) call sync_param( "C_SDRAG", C_SDRAG ) call sync_param( "P_CSDRAG", P_CSDRAG ) call sync_param( "P_SDRAG", P_SDRAG ) call sync_param( "PP_SDRAG", PP_SDRAG ) call sync_param( "ANG_SDRAG", ANG_SDRAG ) call sync_param( "Wc_Jdrag", Wc_Jdrag ) call sync_param( "do_polefix", do_polefix ) call sync_param( "NDASF", NDASF ) call sync_param( "NDA4", NDA4 ) !! call sync_param( "NDA5S", NDA5S ) !! call sync_param( "NDA5K", NDA5K ) !! call sync_param( "NDA5D", NDA5D ) !! call sync_param( "NDAA", NDAA ) !! call sync_param( "NFILTR", NFILTR ) !! call sync_param( "NRAD", NRAD ) !! call sync_param( "Kvflxo", Kvflxo ) !! call sync_param( "Ndisk", Ndisk ) call sync_param( "Nssw", Nssw ) call sync_param( "KCOPY", KCOPY ) call sync_param( "KOCEAN", KOCEAN ) call sync_param( "KRADIA", KRADIA ) call sync_param( "NIsurf", NIsurf ) call sync_param( "UOdrag", UOdrag ) call sync_param( "IRAND", IRAND ) call sync_param( "COUPLED_CHEM", COUPLED_CHEM ) C**** Non-Rundeck parameters C**** Calculate levels for application of SDRAG: LSDRAG,LPSDRAG->LM i.e. C**** all levels above and including P_SDRAG mb (PP_SDRAG near poles) C**** If P is the edge between 2 levels, take the higher level. C**** Also find CSDRAGL, the coefficients of C_Sdrag as a function of L LSDRAG=LM ; LPSDRAG=LM ; LCSDRAG=LM ; CSDRAGL=C_SDRAG DO L=1,LM IF (PEDNL00(L+1)-1d-5.lt.P_SDRAG .and. * PEDNL00(L) +1d-5.gt.P_SDRAG) LSDRAG=L IF (PEDNL00(L+1)-1d-5.lt.PP_SDRAG .and. * PEDNL00(L) +1d-5.gt.PP_SDRAG) LPSDRAG=L IF (PEDNL00(L+1)-1d-5.lt.P_CSDRAG .and. * PEDNL00(L) +1d-5.gt.P_CSDRAG) LCSDRAG=L END DO DO L=LCSDRAG,LSDRAG-1 CSDRAGL(L) = C_SDRAG + max( 0.d0 , (X_SDRAG(1)-C_SDRAG) * * LOG(P_CSDRAG/(PMIDL00(L))) / LOG(P_CSDRAG/P_SDRAG) ) END DO if (AM_I_ROOT()) then WRITE(6,*) "Levels for LSDRAG =",LSDRAG ,"->",LM WRITE(6,*) "Levels for LPSDRAG =",LPSDRAG,"->",LM," near poles" WRITE(6,*) "C_SDRAG coefficients:",CSDRAGL(LS1:LSDRAG-1) end if C**** Determine if FLTRUV is called. QUVfilter = .false. if (DT_XUfilter>0. .or. DT_XVfilter>0. .or. * DT_YUfilter>0. .or. DT_YVfilter>0.) QUVfilter = .true. if (QUVfilter) then if (DT_XUfilter > 0. .and. DT_XUfilter < DT) then DT_XUfilter = DT WRITE(6,*) "DT_XUfilter too small; reset to :",DT_XUfilter end if if (DT_XVfilter > 0. .and. DT_XVfilter < DT) then DT_XVfilter = DT WRITE(6,*) "DT_XVfilter too small; reset to :",DT_XVfilter end if if (DT_YUfilter > 0. .and. DT_YUfilter < DT) then DT_YUfilter = DT WRITE(6,*) "DT_YUfilter too small; reset to :",DT_YUfilter end if if (DT_YVfilter > 0. .and. DT_YVfilter < DT) then DT_YVfilter = DT WRITE(6,*) "DT_YVfilter too small; reset to :",DT_YVfilter end if end if c Warn if polar fixes requested for a model not having a half polar box c if(do_polefix.eq.1 .and. jm.ne.46) then c do_polefix = 0 c write(6,*) 'Polar fixes are currently applicable only to'// c & 'models having a half polar box; no fixes applied' c endif RETURN C**** end subroutine init_Model #ifdef USE_FVCORE SUBROUTINE INPUT (istart,ifile,clock) 2,187 #else SUBROUTINE INPUT (istart,ifile) 2,187 #endif C**** C**** THIS SUBROUTINE SETS THE PARAMETERS IN THE C ARRAY, READS IN THE C**** INITIAL CONDITIONS, AND CALCULATES THE DISTANCE PROJECTION ARRAYS C**** USE FILEMANAGER, only : openunit,closeunit,nameunit USE TIMINGS, only : timing,ntimeacc USE PARAM USE PARSER USE CONSTANT, only : grav,kapa,sday,by3 USE MODEL_COM, only : im,jm,lm,wm,u,v,t,p,q,fearth0,fland * ,focean,flake0,flice,hlake,zatmo,plbot,sig,dsig,sige,kradia * ,bydsig,xlabel,lrunid,nmonav,qcheck,irand,ptop * ,nisurf,nidyn,nday,dt,dtsrc,kdisk,jmon0,jyear0 * ,iyear1,itime,itimei,itimee * ,ls1,psfmpt,pstrat,idacc,jyear,jmon,jday,jdate,jhour * ,aMONTH,jdendofm,jdpery,aMON,aMON0,ioread,irerun * ,ioread_single,irsfic,irsficnt,iowrite_single,ioreadnt * ,irsficno,mdyn,mcnds,mrad,msurf,mdiag,melse,Itime0,Jdate0 * ,Jhour0,rsf_file_name,lm_req * ,pl00,aml00,pednl00,pdsigl00,pmidl00,byaml00 USE SOMTQ_COM, only : tmom,qmom USE GEOM, only : geom_b,imaxj USE RANDOM USE RAD_COM, only : rqt,cloud_rad_forc USE DYNAMICS, only : pk,pmid,pedn USE CLOUDS_COM, only : ttold,qtold,svlhx,rhsav,cldsav #if (defined TRACERS_ON) || (defined TRACERS_OCEAN) USE TRACER_COM,only: MTRACE,NTM,TRNAME #ifdef TRACERS_SPECIAL_Shindell * ,mchem #endif #endif #ifdef TRACERS_AMP USE AERO_CONFIG USE AERO_COAG USE AERO_INIT USE AERO_SETUP USE AERO_SUBS USE AERO_NPF USE AERO_DIAM USE AMP_AEROSOL #endif USE DIAG_COM, only : acc_period,monacc,jreg,titreg,namreg & ,hr_in_day,iwrite,jwrite,itwrite,kdiag,qdiag,qdiag_ratios,oa USE PBLCOM & , only : wsavg,tsavg,qsavg,dclev,usavg,vsavg,tauavg,ustar_pbl & ,egcm,w2gcm,tgvavg,qgavg USE LAKES_COM, only : flake USE GHY_COM, only : fearth USE SOIL_DRV, only: init_gh USE DOMAIN_DECOMP, only : grid, GET, READT_PARALLEL, AM_I_ROOT USE DOMAIN_DECOMP, only : HALO_UPDATE, NORTH, HERE #ifdef USE_FVCORE USE FV_INTERFACE_MOD, only: init_app_clock USE CONSTANT, only : hrday USE ESMF_MOD, only: ESMF_Clock #endif USE ATMDYN, only : init_ATMDYN,CALC_AMPK #ifdef USE_ENT USE ENT_DRV, only : init_module_ent #endif IMPLICIT NONE CHARACTER(*) :: ifile !@var iu_AIC,iu_TOPO,iu_GIC,iu_REG,iu_RSF unit numbers for input files INTEGER iu_AIC,iu_TOPO,iu_GIC,iu_REG,iu_RSF,iu_IFILE !@var num_acc_files number of acc files for diag postprocessing INTEGER I,J,L,K,LID1,LID2,ITYPE,IM1,NOFF,ioerr,num_acc_files !@nlparam HOURI,DATEI,MONTHI,YEARI start of model run !@nlparam TIMEE,HOURE,DATEE,MONTHE,YEARE,IHOURE end of model run !@var IHRI,IHOURE start and end of run in hours (from 1/1/IYEAR1 hr 0) INTEGER :: HOURI=0 , DATEI=1, MONTHI=1, YEARI=-1, IHRI=-1, * TIMEE=-1,HOURE=0 , DATEE=1, MONTHE=1, YEARE=-1, IHOURE=-1, !@nlparam ISTART postprocessing(-1)/start(1-8)/restart(>8) option !@nlparam IRANDI random number seed to perturb init.state (if>0) * ISTART, IRANDI=0 REAL*8 TIJL,CDM,TEMP,X INTEGER Itime1,Itime2,ItimeX,IhrX, LMR !@ egcm_init_max maximum initial vaule of egcm real*8, parameter :: egcm_init_max=0.5d0 LOGICAL :: redoGH = .FALSE.,iniPBL = .FALSE., inilake = .FALSE., & iniSNOW = .FALSE. ! true = restart from "no snow" rsf & ,iniOCEAN = .FALSE. #ifdef USE_ENT & ,iniENT = .FALSE. #endif #ifdef USE_FVCORE type (ESMF_Clock) :: clock integer :: minti,minte character(len=1) :: suffix #endif CHARACTER NLREC*80,filenm*100,RLABEL*132 NAMELIST/INPUTZ/ ISTART,IRANDI * ,IWRITE,JWRITE,ITWRITE,QCHECK,QDIAG,KDIAG,QDIAG_RATIOS * ,IHOURE, TIMEE,HOURE,DATEE,MONTHE,YEARE,IYEAR1 C**** List of parameters that are disregarded at restarts * , HOURI,DATEI,MONTHI,YEARI integer ISTART_kradia, nl_soil c**** Extract domain decomposition info INTEGER :: J_0, J_1, J_0S, J_1S, J_0H, J_1H LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE CCCC INTEGER :: stdin ! used to read 'I' file CALL GET(grid, J_STRT = J_0, J_STOP = J_1, & J_STRT_SKP = J_0S, J_STOP_SKP = J_1S, & J_STRT_HALO= J_0H, J_STOP_HALO = J_1H, & HAVE_SOUTH_POLE = HAVE_SOUTH_POLE, & HAVE_NORTH_POLE = HAVE_NORTH_POLE) C**** C**** Default setting for ISTART : restart from latest save-file (10) C**** ISTART=10 num_acc_files=0 C**** C**** Set dependent vertical resolution variables C**** SIGE(:) = (PLbot(:)-PTOP)/PSFMPT SIG(:) = (sige(1:lm)+sige(2:lm+1))*0.5d0 DSIG(:) = sige(1:lm)-sige(2:lm+1) byDSIG = 1./DSIG C**** CALCULATE SPHERICAL GEOMETRY CALL GEOM_B C**** Calculate default vertical arrays (including rad. eq. layers) LMR=LM+LM_REQ CALL CALC_VERT_AMP(PSFMPT,LMR,PL00,AML00,PDSIGL00,PEDNL00,PMIDL00) BYAML00(:)=1./AML00(:) C**** C**** default settings for prog. variables etc C**** TEMP=250. TSAVG(:,:)=TEMP U(:,:,:)=0. V(:,:,:)=0. T(:,:,:)=TEMP ! will be changed to pot.temp later Q(:,:,:)=3.D-6 P(:,:)=PSFMPT C**** Advection terms for first and second order moments TMOM(:,:,:,:)=0. QMOM(:,:,:,:)=0. C**** Auxiliary clouds arrays RHSAV (:,:,:)=.85d0 CLDSAV(:,:,:)=0. SVLHX (:,:,:)=0. WM (:,:,:)=0. C**** Ocean info saved for ocean heat transport calculations OA = 0. C**** All diagn. are enabled unless KDIAG is changed in the rundeck KDIAG(1:12)=0 KDIAG(13)=9 C**** Set global default timing descriptions C**** Other speciality descriptions can be added/used locally NTIMEACC = 0 CALL SET_TIMER("ATMOS. DYNAM",MDYN) CALL SET_TIMER("CONDENSATION",MCNDS) CALL SET_TIMER(" RADIATION",MRAD) CALL SET_TIMER(" SURFACE",MSURF) CALL SET_TIMER(" DIAGNOSTICS",MDIAG) CALL SET_TIMER(" OTHER",MELSE) #if (defined TRACERS_ON) || (defined TRACERS_OCEAN) CALL SET_TIMER(" TRACERS",MTRACE) #endif #ifdef TRACERS_SPECIAL_Shindell CALL SET_TIMER(" CHEMISTRY",MCHEM) #endif C**** C**** Set some documentary parameters in the database C**** call set_param("IM",IM) call set_param("JM",JM) call set_param("LM",LM) call set_param("LS1",LS1) call set_param("PLBOT",Plbot,LM+1) #if (defined TRACERS_ON) || (defined TRACERS_OCEAN) call set_param("NTM",NTM) call set_param("TRNAME",TRNAME,ntm) #endif C**** C**** Print Header and Label (2 lines) from rundeck C**** call openunit(trim(ifile),iu_IFILE,.false.,.true.) if (AM_I_ROOT()) WRITE (6,'(A,40X,A/)') '0','GISS CLIMATE MODEL' READ(iu_IFILE,'(A80)') XLABEL(1:80),NLREC NOFF=0 IF (XLABEL(73:80).EQ.' ') NOFF=8 ! for 72-column rundecks XLABEL(81-NOFF:132)=NLREC(1:52+NOFF) if (AM_I_ROOT()) WRITE (6,'(A,A/)') '0',XLABEL RLABEL = XLABEL !@var RLABEL rundeck-label IF(AM_I_ROOT()) THEN C**** C**** Print preprocessing options (if any are defined) C**** #if (defined TRACERS_ON) || (defined TRACERS_OCEAN) write(6,*) 'This program includes tracer code' #endif #ifdef TRACERS_WATER write(6,*) '...and water tracer code' #ifndef TRACERS_ON call stop_model( &' Water tracers need TRACERS_ON as well as TRACERS_WATER',255) #endif #endif #ifdef TRACERS_OCEAN write(6,*) '...and ocean tracer code' #endif #ifdef TRACERS_SPECIAL_O18 write(6,*) '...and water isotope code' #ifndef TRACERS_WATER call stop_model('Water isotope tracers need TRACERS_WATER '// *'as well as TRACERS_SPECIAL_O18',255) #endif #endif #ifdef TRACERS_SPECIAL_Lerner write(6,*) '...and Jean/David tracers and chemistry' #endif #ifdef TRACERS_GASEXCH_Natassa write(6,*) ' ' write(6,*) '...and Natassa Romanou air-sea GAS EXCHANGE' #ifdef TRACERS_GASEXCH_CFC_Natassa write(6,*) '****CFC flux across air/sea interface****' #endif #ifdef TRACERS_GASEXCH_CO2_Natassa write(6,*) '****CO2 flux across air/sea interface****' #endif #endif #ifdef TRACERS_SPECIAL_Shindell write(6,*) '...and Drew Shindell tracers and chemistry' #endif #ifdef TRACERS_AEROSOLS_Koch write(6,*) '...and Dorothy Koch aerosols' #endif #ifdef TRACERS_DRYDEP write(6,*) '...and tracer dry deposition' #endif #ifdef EDGAR_HYDE_SOURCES write(6,*) '...and EDGAR HYDE sources instead of GISS' #endif #ifdef SHINDELL_STRAT_CHEM write(6,*) '...and Drew Shindell stratospheric chemistry' #endif #ifdef SHINDELL_STRAT_EXTRA write(6,*) '...and Drew Shindell extra strat tracers' #endif #ifdef regional_Ox_tracers write(6,*) '...and regional Ox tracers' #endif #ifdef INTERACTIVE_WETLANDS_CH4 write(6,*) '...and interactive CH4 wetlands emissions' #endif #ifdef NUDGE_ON write(6,*) '...and nudging of meteorology' #endif ENDIF ! AM_I_ROOT() C**** C**** Read parameters from the rundeck to database and namelist C**** call parse_params(iu_IFILE) READ (iu_IFILE,NML=INPUTZ,ERR=900) call closeunit(iu_IFILE) C**** Get those parameters which are needed in this subroutine if(is_set_param("DTsrc")) call get_param( "DTsrc", DTsrc ) if(is_set_param("DT")) call get_param( "DT", DT ) if(is_set_param("NIsurf")) call get_param( "NIsurf", NIsurf ) ! if(is_set_param("IRAND")) call get_param( "IRAND", IRAND ) if(is_set_param("NMONAV")) call get_param( "NMONAV", NMONAV ) if(is_set_param("Kradia")) call get_param( "Kradia", Kradia ) C*********************************************************************** C**** **** C**** Post-process one or more ACC-files : ISTART < 1 **** C**** **** C*********************************************************************** if (istart.le.0) then call reset_diag(1) monacc = 0 do call nextarg(filenm, 0) if ( filenm == "" ) exit ! end of args call openunit(filenm,iu_AIC,.true.,.true.) call io_rsf(iu_AIC,itime,ioread_single,ioerr) call closeunit(iu_AIC) num_acc_files = num_acc_files + 1 end do GO TO 500 end if if (istart.ge.9 .or. Kradia.gt.0) go to 400 C*********************************************************************** C**** **** C**** INITIAL STARTS - ISTART: 1 to 8 **** C**** **** C**** Current settings: 1 - from defaults **** C**** 2 - from observed data **** C**** 3 - so far unused **** C**** 4 - from coupled model M-file - reset ocn **** C**** 5 - tracer run from M-file w/o tracers **** C**** 6 - pred.ocn run from M-file w/o ocn data **** C**** 7 - from mod. II' M-file - reset snow/ocn **** C**** 8 - from current model M-file - no resets **** C**** **** C*********************************************************************** C**** get unit for atmospheric initial conditions if needed IF (ISTART.gt.1) call openunit("AIC",iu_AIC,.true.,.true.) C**** C**** Set quantities that are derived from the namelist parameters C**** !@var NDAY=(1 day)/DTsrc : even integer; adjust DTsrc later if necessary NDAY = 2*NINT(.5*SDAY/DTsrc) C**** Get Start Time; at least YearI HAS to be specified in the rundeck IF (YearI.lt.0) then IF (AM_I_ROOT()) * WRITE(6,*) 'Please choose a proper start year yearI, not',yearI call stop_model('INPUT: yearI not provided',255) END IF IF (Iyear1.lt.0) Iyear1 = yearI IhrI = HourI + + HR_IN_DAY*(dateI-1 + JDendofM(monthI-1) + JDperY*(yearI-Iyear1)) ITimeI = IhrI*NDAY/HR_IN_DAY ! internal clock counts DTsrc-steps Itime=ItimeI IF (IhrI.lt.0) then IF (AM_I_ROOT()) * WRITE(6,*) 'Improper start time OR Iyear1=',Iyear1,' > yearI;', * ' yearI,monthI,dateI,hourI=',yearI,monthI,dateI,hourI call stop_model( & 'INPUT: Improper start date or base year Iyear1',255) END IF C**** Check the vertical layering defined in RES_ (is sige(ls1)=0 ?) IF (SIGE(LS1).ne.0.) then if (AM_I_ROOT()) * write(6,*) 'bad vertical layering: ls1,sige(ls1)', & ls1,sige(ls1) call stop_model('INPUT: ls1 incorrectly set in RES_',255) END IF C**** C**** Get Ground conditions from a separate file - ISTART=1,2 C**** IF (ISTART.LE.2) THEN C**** Set flag to initialise pbl and snow variables iniPBL=.TRUE. iniSNOW = .TRUE. ! extract snow data from first soil layer iniOCEAN = .TRUE. ! read in ocean ic #ifdef USE_ENT iniENT = .TRUE. #endif if (istart.eq.1) redogh=.true. C**** Read in ground initial conditions call openunit("GIC",iu_GIC,.true.,.true.) ioerr=-1 read(iu_GIC) ! ignore first line (ocean ic done in init_OCEAN) call io_seaice (iu_GIC,ioreadnt,ioerr) call io_earth (iu_GIC,ioreadnt,ioerr) call io_soils (iu_GIC,ioreadnt,ioerr) call io_landice(iu_GIC,ioreadnt,ioerr) if (ioerr.eq.1) then IF (AM_I_ROOT()) * WRITE(6,*) "I/O ERROR IN GIC FILE: KUNIT=",iu_GIC call stop_model("INPUT: GIC READ IN ERROR",255) end if call closeunit (iu_GIC) END IF C**** C**** Get primary Atmospheric data from NMC tapes - ISTART=2 C**** IF (ISTART.EQ.2) THEN C**** Use title of first record to get the date and make sure ??? C**** it is consistent with IHRI (at least equal mod 8760) ??? C**** not yet implemented but could easily be done ??? XLABEL(1:80)='Observed atmospheric data from NMC tape' Csoon READ (iu_AIC) XLABEL(1:80) CALL READT_PARALLEL(grid,iu_AIC,NAMEUNIT(iu_AIC),0,P,1) ! Psurf DO J=J_0,J_1 DO I=1,IM P(I,J)=P(I,J)-PTOP ! Psurf -> P END DO END DO DO L=1,LM CALL READT_PARALLEL(grid,iu_AIC,NAMEUNIT(iu_AIC),0,U(:,:,L),1) ! U END DO DO L=1,LM CALL READT_PARALLEL(grid,iu_AIC,NAMEUNIT(iu_AIC),0,V(:,:,L),1) ! V END DO DO L=1,LM CALL READT_PARALLEL(grid,iu_AIC,NAMEUNIT(iu_AIC),0,T(:,:,L),1) ! Temperature END DO DO L=1,LM ! alternatively, only read in L=1,LS1 ; skip rest CALL READT_PARALLEL(grid,iu_AIC,NAMEUNIT(iu_AIC),0,Q(:,:,L),1) ! Q END DO CALL READT_PARALLEL(grid,iu_AIC,NAMEUNIT(iu_AIC),0,TSAVG,1) ! Tsurf END IF C**** C**** Derive other data from primary data if necessary - ISTART=1,2 C**** currently IF (ISTART.LE.2) THEN If (HAVE_SOUTH_POLE) THEN WSAVG(1:im,1)=SQRT(U(1,2,1)*U(1,2,1)+V(1,2,1)*V(1,2,1)) USAVG(1:im,1)=U(1,2,1) VSAVG(1:im,1)=V(1,2,1) End If If (HAVE_NORTH_POLE) THEN WSAVG(1:im,JM)=SQRT(U(1,JM,1)*U(1,JM,1)+V(1,JM,1)*V(1,JM,1)) USAVG(1:im,JM)=U(1,JM,1) VSAVG(1:im,JM)=V(1,JM,1) End If CALL HALO_UPDATE(grid, U, FROM=NORTH) CALL HALO_UPDATE(grid, V, FROM=NORTH) DO J=J_0S,J_1S IM1=IM DO I=1,IM WSAVG(I,J)=.25*SQRT( * (U(IM1,J,1)+U(I,J,1)+U(IM1,J+1,1)+U(I,J+1,1))**2 * +(V(IM1,J,1)+V(I,J,1)+V(IM1,J+1,1)+V(I,J+1,1))**2) USAVG(I,J)=.25*(U(IM1,J,1)+U(I,J,1)+U(IM1,J+1,1)+U(I,J+1,1)) VSAVG(I,J)=.25*(V(IM1,J,1)+V(I,J,1)+V(IM1,J+1,1)+V(I,J+1,1)) IM1=I END DO END DO CDM=.001d0 CALL CALC_AMPK(LM) DO J=J_0,J_1 DO I=1,IM C**** SET SURFACE MOMENTUM TRANSFER TAU0 TAUAVG(I,J)=CDM*WSAVG(I,J)**2 C**** SET LAYER THROUGH WHICH DRY CONVECTION MIXES TO 1 DCLEV(I,J)=1. C**** SET SURFACE SPECIFIC HUMIDITY FROM FIRST LAYER HUMIDITY QSAVG(I,J)=Q(I,J,1) QGAVG(I,J)=Q(I,J,1) TGVAVG(I,J)=T(I,J,1) C**** SET RADIATION EQUILIBRIUM TEMPERATURES FROM LAYER LM TEMPERATURE DO K=1,LM_REQ RQT(K,I,J)=T(I,J,LM) END DO C**** REPLACE TEMPERATURE BY POTENTIAL TEMPERATURE DO L=1,LM T(I,J,L)=T(I,J,L)/PK(L,I,J) TTOLD(L,I,J)=T(I,J,L) QTOLD(L,I,J)=Q(I,J,L) END DO C**** initialize egcm to be used in ATURB.f DO L=1,LM egcm(l,i,j)=egcm_init_max/(float(l)**2) w2gcm(l,i,j)=egcm(l,i,j)*2.*by3 END DO END DO END DO C**** Initialize surface friction velocity DO ITYPE=1,4 DO J=J_0,J_1 DO I=1,IM USTAR_pbl(I,J,ITYPE)=WSAVG(I,J)*SQRT(CDM) END DO END DO END DO C**** INITIALIZE VERTICAL SLOPES OF T,Q call tq_zmom_init(t,q,PMID,PEDN) END IF C**** C**** I.C from possibly older/incomplete MODEL OUTPUT, ISTART=3-8 C**** SELECT CASE (ISTART) CASE (3) ! just general hints - not to be used as is C**** Read what's there and substitute rest as needed (as above) C**** To be implemented as needed. Sometimes it is safer to C**** combine the ground layers into 2 layers (top 10cm and rest) and C**** set redoGH to .true. (after major changes in the GH code or C**** after changing to a new horizontal grid) C redoGH=.TRUE. C**** Set flag to initialise pbl/snow variables if obsolete or missing C iniPBL=.TRUE. ; iniSNOW = .TRUE. go to 890 ! not implemented; stop model C**** C**** I.C FROM FULL MODEL RESTART FILE (but re-initialise ocean) C**** CASE (4) call io_rsf(iu_AIC,IhrX,irsficno,ioerr) if (ioerr.eq.1) goto 800 iniOCEAN = .TRUE. ! read in ocean ic C**** C**** I.C FROM FULL MODEL RESTART FILE (but no tracers) C**** CASE (5) ! this model's rsf file, no tracers call io_rsf(iu_AIC,IhrX,irsficnt,ioerr) if (ioerr.eq.1) goto 800 C**** C**** I.C FROM RESTART FILE that may not match land-ocean mask ISTART=6 C**** CASE (6) ! converted model II' (B399) format (no snow) call io_rsf(iu_AIC,IhrX,irsficno,ioerr) if (ioerr.eq.1) goto 800 iniSNOW = .TRUE. ! extract snow data from first soil layer inipbl = .TRUE. ! initialise pbl profiles iniOCEAN = .TRUE. ! read in ocean ic redoGH=.TRUE. C**** C**** I.C FROM RESTART FILE WITH almost COMPLETE DATA ISTART=7 C**** CASE (7) ! converted model II' (B399) format (no snow) call io_rsf(iu_AIC,IhrX,irsfic,ioerr) if (ioerr.eq.1) goto 800 iniSNOW = .TRUE. ! extract snow data from first soil layer iniOCEAN = .TRUE. ! read in ocean ic C**** C**** Data from current type of RESTART FILE ISTART=8 C**** CASE (8) ! no need to read SRHR,TRHR,FSF,TSFREZ,diag.arrays call io_rsf(iu_AIC,IhrX,irsfic,ioerr) iniSNOW = .TRUE. ! extract snow data from first soil layer iniPBL=.TRUE. if (ioerr.eq.1) goto 800 END SELECT C**** Check consistency of starting time IF (ISTART.ge.3) THEN IF( (MOD(IHRI-IHRX,8760).ne.0) ) THEN WRITE (6,*) ' Difference in hours between ', * 'Starting date and Data date:',MOD(IHRI-IHRX,8760) WRITE (6,*) 'Please change HOURI,DATEI,MONTHI' call stop_model('INPUT: start date inconsistent with data',255) ENDIF END IF C**** Set flag to initialise lake variables if they are not in I.C. IF (ISTART.lt.8) inilake=.TRUE. CALL CALC_AMPK(LM) C**** !**** IRANDI seed for random perturbation of initial conditions (if/=0): C**** tropospheric temperatures changed by at most 1 degree C IF (IRANDI.NE.0) THEN CALL RINIT (IRANDI) DO L=1,LS1-1 call burn_random(n=im*(j_0-1)) DO J=J_0,J_1 DO I=1,IM TIJL=T(I,J,L)*PK(L,I,J)-1.+2*RANDU(X) T(I,J,L)=TIJL/PK(L,I,J) END DO END DO END DO IF (AM_I_ROOT()) * WRITE(6,*) 'Initial conditions were perturbed !!',IRANDI END IF C**** Close "AIC" here if it was opened IF (ISTART.gt.1) call closeunit(iu_AIC) IF (AM_I_ROOT()) * WRITE(6,'(A,i3,1x,a4,i5,a3,i3,3x,a,i2/" ",a)') * '0Model started on',datei,aMONTH(monthi),yeari,' Hr',houri, * 'ISTART =',ISTART,XLABEL(1:80) ! report input file label XLABEL = RLABEL ! switch to rundeck label GO TO 600 C*********************************************************************** C**** **** C**** RESTARTS: ISTART > 8 **** C**** **** C**** Current settings: 9 - from own model M-file **** C**** 10 - from later of fort.1 or fort.2 **** C**** 11 - from fort.1 **** C**** 12 - from fort.2 **** C**** 13 & up - from earlier of fort.1 or fort.2 **** C**** **** C*********************************************************************** 400 SELECT CASE (ISTART) C**** C**** DATA FROM end-of-month RESTART FILE ISTART=9 C**** mainly used for REPEATS and delayed EXTENSIONS CASE (1:9) ! diag.arrays are not read in call openunit("AIC",iu_AIC,.true.,.true.) if(istart.eq.9) call io_rsf(iu_AIC,Itime,irerun,ioerr) #ifdef USE_FVCORE call system('cp AICfv fv_internal_restart.dat') call system('cp AICdfv tendencies_checkpoint') #endif if(istart.le.8) then ! initial start of rad.forcing run call io_label(iu_AIC,Itime,ItimeX,irerun,ioerr) if (Kradia.gt.0) call io_rad (iu_AIC,irsfic,ioerr) end if call closeunit(iu_AIC) if (ioerr.eq.1) goto 800 WRITE (6,'(A,I2,A,I11,A,A/)') '0Model restarted; ISTART=', * ISTART,', TIME=',Itime,' ',XLABEL(1:80) ! sho input file label XLABEL = RLABEL ! switch to rundeck label CALL CALC_AMPK(LM) C**** !**** IRANDI seed for random perturbation of current state (if/=0) C**** tropospheric temperatures are changed by at most 1 degree C IF (IRANDI.ne.0 .and. Kradia.le.0) THEN CALL RINIT (IRANDI) DO L=1,LS1-1 DO J=J_0,J_1 DO I=1,IM TIJL=T(I,J,L)*PK(L,I,J)-1.+2*RANDU(X) T(I,J,L)=TIJL/PK(L,I,J) END DO END DO END DO IF (AM_I_ROOT()) * WRITE(6,*) 'Current temperatures were perturbed !!',IRANDI END IF TIMING = 0 GO TO 500 C**** C**** RESTART ON DATA SETS 1 OR 2, ISTART=10 or more C**** C**** CHOOSE DATA SET TO RESTART ON CASE (10,13:) Itime1=-1 call openunit(rsf_file_name(1),iu_RSF,.true.,.true.) READ (iu_RSF,ERR=410) Itime1 call closeunit(iu_RSF) 410 continue !REWIND 1 Itime2=-1 call openunit(rsf_file_name(2),iu_RSF,.true.,.true.) READ (iu_RSF,ERR=420) Itime2 call closeunit(iu_RSF) 420 continue !REWIND 2 IF (Itime1+Itime2.LE.-2.) GO TO 850 KDISK=1 IF (Itime2.GT.Itime1) KDISK=2 IF (ISTART.GE.13) KDISK=3-KDISK CASE (11,12) KDISK=ISTART-10 END SELECT 430 continue call openunit(rsf_file_name(KDISK),iu_RSF,.true.,.true.) #ifdef USE_FVCORE write(suffix,'(i1)') kdisk call system('cp fv.'// suffix // ' fv_internal_restart.dat') call system('cp dfv.'// suffix // ' tendencies_checkpoint') #endif CALL HERE(__FILE__//'::io_rsf',__LINE__ + 10000*KDISK) call io_rsf(iu_RSF,Itime,ioread,ioerr) call closeunit(iu_RSF) if (ioerr.eq.1) then if (istart.gt.10) go to 850 ! no 2nd chance if istart/=10 KDISK=3-KDISK ! try the earlier restart file WRITE (6,'(A,I1,A,I1)') * ' Read Error on fort.',3-kdisk,' trying fort.',kdisk ISTART=110 go to 430 end if if (AM_I_ROOT()) * WRITE (6,'(A,I2,A,I11,A,A/)') '0RESTART DISK READ, UNIT', * KDISK,', Time=',Itime,' ',XLABEL(1:80) C**** Switch KDISK if the other file is (or may be) bad (istart>10) C**** so both files will be fine after the next write execution IF (istart.gt.10) KDISK=3-KDISK C**** Keep KDISK after reading from the later restart file, so that C**** the same file is overwritten first; in case of trouble, C**** the earlier restart file will still be available 500 CONTINUE C**** Get parameters we just read from rsf file. Only those C**** parameters which we need in "INPUT" should be extracted here. if(is_set_param("DTsrc")) call get_param( "DTsrc", DTsrc ) if(is_set_param("DT")) call get_param( "DT", DT ) if(is_set_param("NMONAV")) call get_param( "NMONAV", NMONAV ) if(is_set_param("Kradia")) call get_param( "Kradia", Kradia ) C*********************************************************************** C**** ***** C**** INITIAL- AND RESTARTS: Final Initialization steps ***** C**** ***** C*********************************************************************** 600 CONTINUE C**** initialize Lrunid (length of the identifying part of XLABEL) C**** lid1 = INDEX(XLABEL,'(') -1 if (lid1.lt.1) lid1=17 lid2 = INDEX(XLABEL,' ') -1 if (lid2.lt.1) lid2=17 LRUNID = min(lid1,lid2) IF (LRUNID.gt.16) call stop_model * ('INPUT: Rundeck name too long. Shorten to 16 char or less' * ,255) C**** Update ItimeE only if YearE or IhourE is specified in the rundeck C**** if(timee.lt.0) timee=houre*nday/HR_IN_DAY IF(yearE.ge.0) ItimeE = (( (yearE-iyear1)*JDperY + * JDendofM(monthE-1)+dateE-1 )*HR_IN_DAY )*NDAY/HR_IN_DAY + TIMEE C**** Alternate (old) way of specifying end time if(IHOURE.gt.0) ItimeE=IHOURE*NDAY/HR_IN_DAY C**** Check consistency of DTsrc (with NDAY) and dt (with NIdyn) if (is_set_param("DTsrc") .and. nint(sday/DTsrc).ne.NDAY) then if (AM_I_ROOT()) * write(6,*) 'DTsrc=',DTsrc,' has to stay at/be set to',SDAY/NDAY call stop_model('INPUT: DTsrc inappropriately set',255) end if DTsrc = SDAY/NDAY call set_param( "DTsrc", DTsrc, 'o' ) ! copy DTsrc into DB NIdyn=nint(dtsrc/dt) #ifndef USE_FVCORE C**** NIdyn=dtsrc/dt(dyn) has to be a multiple of 2 C**** if(istart>0) then NIdyn = 2*nint(.5*dtsrc/dt) if (is_set_param("DT") .and. nint(DTsrc/dt).ne.NIdyn) then if (AM_I_ROOT()) * write(6,*) 'DT=',DT,' has to be changed to',DTsrc/NIdyn call stop_model('INPUT: DT inappropriately set',255) end if end if #else C**** need a clock to satisfy ESMF interfaces call getdte(itimei,nday,iyear1,YEARI,MONTHI,jday,DATEI,HOURI,amon) MINTI = nint(mod( mod(Itimei*hrday/Nday,hrday) * 60d0, 60d0)) call getdte(itimee,nday,iyear1,YEARE,MONTHE,jday,DATEE,HOURE,amon) MINTE = nint(mod( mod(Itimee*hrday/Nday,hrday) * 60d0, 60d0)) clock = init_app_clock( (/ YEARI, MONTHI, DATEI, HOURI, MINTI,0/), & (/ YEARE, MONTHE, DATEE, HOURE, MINTE,0/), & interval = int(dt) ) #endif DT = DTsrc/NIdyn call set_param( "DT", DT, 'o' ) ! copy DT into DB C**** NMONAV has to be 1(default),2,3,4,6,12, i.e. a factor of 12 if (NMONAV.lt.1 .or. MOD(12,NMONAV).ne.0) then write (6,*) 'NMONAV has to be 1,2,3,4,6 or 12, not',NMONAV call stop_model('INPUT: nmonav inappropriately set',255) end if if (AM_I_ROOT()) * write (6,*) 'Diag. acc. period:',NMONAV,' month(s)' C**** Updating Parameters: If any of them changed beyond this line C**** use set_param(.., .., 'o') to update them in the database (DB) C**** Get the rest of parameters from DB or put defaults to DB call init_Model C**** Set julian date information call getdte(Itime,Nday,Iyear1,Jyear,Jmon,Jday,Jdate,Jhour,amon) call getdte(Itime0,Nday,iyear1,Jyear0,Jmon0,J,Jdate0,Jhour0,amon0) C**** C**** READ IN TIME-INDEPENDENT ARRAYS C**** if (Kradia.le.0) then ! full model CALL CALC_AMPK(LM) C**** READ SPECIAL REGIONS FROM UNIT 29 call openunit("REG",iu_REG,.true.,.true.) READ(iu_REG) TITREG,JREG,NAMREG IF (AM_I_ROOT()) & WRITE(6,*) ' read REGIONS from unit ',iu_REG,': ',TITREG call closeunit(iu_REG) end if ! full model: Kradia le 0 #if (defined TRACERS_ON) || (defined TRACERS_OCEAN) C**** Initialise tracer parameters and diagnostics call init_tracer #endif C**** READ IN LANDMASKS AND TOPOGRAPHIC DATA C**** Note that FLAKE0 is read in only to provide initial values C**** Actual array is set from restart file. call openunit("TOPO",iu_TOPO,.true.,.true.) CALL READT_PARALLEL(grid,iu_TOPO,NAMEUNIT(iu_TOPO),0,FOCEAN,1) ! Ocean fraction CALL READT_PARALLEL(grid,iu_TOPO,NAMEUNIT(iu_TOPO),0,FLAKE0,1) ! Orig. Lake fraction CALL READT_PARALLEL(grid,iu_TOPO,NAMEUNIT(iu_TOPO),0,FEARTH0,1) ! Earth frac. (no LI) CALL READT_PARALLEL(grid,iu_TOPO,NAMEUNIT(iu_TOPO),0,FLICE ,1) ! Land ice fraction CALL READT_PARALLEL(grid,iu_TOPO,NAMEUNIT(iu_TOPO),0,ZATMO ,1) ! Topography CALL READT_PARALLEL(grid,iu_TOPO,NAMEUNIT(iu_TOPO),0,HLAKE ,2) ! Lake Depths ZATMO(:,J_0:J_1) = ZATMO(:,J_0:J_1)*GRAV ! Geopotential call closeunit(iu_TOPO) C**** Initialise some modules before finalising Land/Ocean/Lake/LI mask C**** Initialize ice CALL init_ice(iniOCEAN,istart) C**** Initialize lake variables (including river directions) CALL init_LAKES(inilake,istart) C**** Initialize ocean variables C**** KOCEAN = 1 => ocean heat transports/max. mixed layer depths C**** KOCEAN = 0 => RSI/MSI factor CALL init_OCEAN(iniOCEAN,istart) C**** Initialize ice dynamics code (if required) CALL init_icedyn(iniOCEAN) C**** Initialize land ice (must come after oceans) CALL init_LI(istart) C**** Make sure that constraints are satisfied by defining FLAND/FEARTH C**** as residual terms. (deals with SP=>DP problem) DO J=J_0,J_1 DO I=1,IMAXJ(J) IF (FOCEAN(I,J).gt.0) THEN FLAND(I,J)=1.-FOCEAN(I,J) ! Land fraction IF (FLAKE(I,J).gt.0) THEN IF (AM_I_ROOT()) WRITE(6,*) * "Ocean and lake cannot co-exist in same grid box" * ,i,j,FOCEAN(I,J),FLAKE(I,J) FLAKE(I,J)=0 END IF ELSEIF (FLAKE(I,J).gt.0) THEN FLAND(I,J)=1.-FLAKE(I,J) ELSE FLAND(I,J)=1. END IF C**** Ensure that no round off error effects land with ice and earth IF (FLICE(I,J)-FLAND(I,J).gt.-1d-4 .and. FLICE(I,J).gt.0) THEN FLICE(I,J)=FLAND(I,J) FEARTH(I,J)=0. ELSE FEARTH(I,J)=FLAND(I,J)-FLICE(I,J) ! Earth fraction END IF END DO END DO If (HAVE_SOUTH_POLE) Then FLAND(2:IM,1)=FLAND(1,1) FEARTH(2:IM,1)=FEARTH(1,1) FLICE(2:IM,1)=FLICE(1,1) End If If (HAVE_NORTH_POLE) Then FLAND(2:IM,JM)=FLAND(1,JM) FEARTH(2:IM,JM)=FEARTH(1,JM) FLICE(2:IM,JM)=FLICE(1,JM) End If C**** C**** INITIALIZE GROUND HYDROLOGY ARRAYS (INCL. VEGETATION) C**** Recompute Ground hydrology data if redoGH (new soils data) C**** !!! hack: make sure that ISTART_kradia==0 if Kradia>0 !!! do we need it ? I.A. ISTART_kradia = ISTART if ( Kradia.gt.0 ) ISTART_kradia = 0 CALL init_GH(DTsrc/NIsurf,redoGH,iniSNOW,ISTART_kradia, nl_soil) #ifdef USE_ENT CALL init_module_ent(iniENT, Jday, Jyear, FOCEAN) !!! FEARTH) #endif if (Kradia.gt.0) then ! radiative forcing run !CALL init_GH(DTsrc/NIsurf,redoGH,iniSNOW,0) !CALL init_module_ent(iniENT,grid,jday,dxyp) CALL init_RAD(istart) if(istart.lt.0) CALL init_DIAG(0,num_acc_files) !post-processing if (AM_I_ROOT()) Then WRITE (6,INPUTZ) call print_param( 6 ) WRITE (6,'(A14,4I4)') "IM,JM,LM,LS1=",IM,JM,LM,LS1 WRITE (6,*) "PLbot=",PLbot end if if(istart.lt.0) & CALL stop_model ('Terminated normally, istart<0',13) return end if ! Kradia>0; radiative forcing run ! CALL init_GH(DTsrc/NIsurf,redoGH,iniSNOW,ISTART) ! CALL init_module_ent(iniENT,grid,jday,dxyp) C**** Initialize pbl (and read in file containing roughness length data) if(istart.gt.0) CALL init_pbl(iniPBL) C**** C**** Initialize the use of gravity wave drag diagnostics C**** CALL init_GWDRAG C C**** Initialize nudging #ifdef NUDGE_ON CALL NUDGE_INIT #endif #ifdef TRACERS_AMP CALL SETUP_CONFIG CALL SETUP_SPECIES_MAPS CALL SETUP_DP0 CALL SETUP_AERO_MASS_MAP CALL SETUP_COAG_TENSORS CALL SETUP_DP0 CALL SETUP_KIJ CALL SETUP_EMIS CALL SETUP_KCI CALL SETUP_NPFMASS CALL SETUP_DIAM #endif C**** if(istart.gt.0) CALL RINIT (IRAND) CALL FFT0 (IM) CALL init_CLD CALL init_DIAG(istart,num_acc_files) ! initialize for accumulation CALL UPDTYPE if(istart.gt.0) CALL init_QUS(grid,im,jm,lm) if(istart.gt.0) CALL init_ATMDYN CALL init_RAD(istart) if (AM_I_ROOT()) then WRITE (6,INPUTZ) call print_param( 6 ) WRITE (6,'(A7,12I6)') "IDACC=",(IDACC(I),I=1,12) WRITE (6,'(A14,4I4)') "IM,JM,LM,LS1=",IM,JM,LM,LS1 WRITE (6,*) "PLbot=",PLbot end if C**** RETURN C**** C**** TERMINATE BECAUSE OF IMPROPER PICK-UP C**** 800 WRITE (6,'(A,I4/" ",A)') * '0ERROR ENCOUNTERED READING AIC ISTART=', ISTART,XLABEL(1:80) call stop_model('INPUT: READ ERROR FOR AIC',255) 830 WRITE(6,*) 'READ ERROR FOR GIC' call stop_model('INPUT: READ ERROR FOR GIC',255) 850 WRITE (6,'(A)') * '0ERRORS ON BOTH RESTART DATA SETS. TERMINATE THIS JOB' call stop_model('INPUT: ERRORS ON BOTH RESTART FILES',255) 890 WRITE (6,'(A,I5)') '0INCORRECT VALUE OF ISTART',ISTART call stop_model('INPUT: ISTART-SPECIFICATION INVALID',255) 900 write (6,*) 'Error in NAMELIST parameters' call stop_model('Error in NAMELIST parameters',255) END SUBROUTINE INPUT SUBROUTINE DAILY(end_of_day) 3,14 !@sum DAILY performs daily tasks at end-of-day and maybe at (re)starts !@auth Original Development Team !@ver 1.0 !@calls constant:orbit, calc_ampk, getdte USE MODEL_COM, only : im,jm,lm,ls1,ptop,psf,p,q * ,itime,itimei,iyear1,nday,jdpery,jdendofm * ,jyear,jmon,jday,jdate,jhour,aMON,aMONTH,ftype USE GEOM, only : areag,dxyp,imaxj USE DYNAMICS, only : byAM USE RADPAR, only : ghgam,ghgyr2,ghgyr1 USE RAD_COM, only : RSDIST,COSD,SIND, dh2o,H2ObyCH4,ghg_yr, * omegt,obliq,eccn #ifdef TRACERS_WATER USE TRACER_COM, only: trm,tr_wd_type,nwater,tr_H2ObyCH4,itime_tr0 * ,ntm #endif USE DIAG_COM, only : aj=>aj_loc,j_h2och4 USE DOMAIN_DECOMP, only : grid, GET, GLOBALSUM, AM_I_ROOT USE ATMDYN, only : CALC_AMPK IMPLICIT NONE REAL*8 DELTAP,PBAR,SMASS,LAM,xCH4,EDPY,VEDAY REAL*8 :: SPRESS(grid%J_STRT_HALO:grid%J_STOP_HALO) INTEGER i,j,l,iy LOGICAL, INTENT(IN) :: end_of_day #ifdef TRACERS_WATER INTEGER n #endif c**** Extract domain decomposition info INTEGER :: J_0, J_1, J_0S, J_1S LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE CALL GET(grid, J_STRT = J_0, J_STOP = J_1, & J_STRT_SKP = J_0S, J_STOP_SKP = J_1S, & HAVE_SOUTH_POLE = HAVE_SOUTH_POLE, & HAVE_NORTH_POLE = HAVE_NORTH_POLE) C**** Tasks to be done at end of day and at each start or restart C**** C**** CALCULATE THE DAILY CALENDAR C**** call getdte(Itime,Nday,iyear1,Jyear,Jmon,Jday,Jdate,Jhour,amon) C**** CALCULATE SOLAR ANGLES AND ORBIT POSITION C**** This is for noon (GMT) for new day. C**** The orbital calculation will need to vary depending on the kind C**** of calendar adopted (i.e. a generic 365 day year, or a transient C**** calendar including leap years etc.). For transient calendars the C**** JDAY passed to orbit needs to be adjusted to represent the number C**** of days from Jan 1 2000AD. c EDPY=365.2425d0, VEDAY=79.3125d0 ! YR 2000AD c JDAY => JDAY + 365 * (JYEAR-2000) + appropriate number of leaps C**** Default calculation (no leap, VE=Mar 21 hr 0) EDPY=365d0 ; VEDAY=79d0 ! Generic year CALL ORBIT (OBLIQ,ECCN,OMEGT,VEDAY,EDPY,REAL(JDAY,KIND=8)-.5 * ,RSDIST,SIND,COSD,LAM) IF (.not.(end_of_day.or.itime.eq.itimei)) RETURN C**** Tasks to be done at end of day and at initial starts only C**** C**** THE GLOBAL MEAN PRESSURE IS KEPT CONSTANT AT PSF MILLIBARS C**** C**** CALCULATE THE CURRENT GLOBAL MEAN PRESSURE SMASS=0. DO J=J_0,J_1 SPRESS(J)=0. DO I=1,IM SPRESS(J)=SPRESS(J)+P(I,J) END DO SPRESS(J) = SPRESS(J) * DXYP(J) END DO CALL GLOBALSUM(grid, SPRESS, SMASS, ALL=.TRUE.) PBAR=SMASS/AREAG+PTOP C**** CORRECT PRESSURE FIELD FOR ANY LOSS OF MASS BY TRUNCATION ERROR C**** except if it was just done (restart from itime=itimei) DELTAP=PSF-PBAR if(itime.eq.itimei .and. abs(deltap).lt.1.d-10) return P=P+DELTAP CALL CALC_AMPK(LS1-1) if (AM_I_ROOT()) then IF (ABS(DELTAP).gt.1d-6) * WRITE (6,'(A25,F10.6/)') '0PRESSURE ADDED IN GMP IS',DELTAP end if IF (.not.end_of_day) RETURN C**** Tasks to be done at end of day only if (H2ObyCH4.gt.0) then C**** Add obs. H2O generated by CH4(*H2ObyCH4) using a 2 year lag iy = jyear - 2 - ghgyr1 + 1 if (ghg_yr.gt.0) iy = ghg_yr - 2 - ghgyr1 + 1 if (iy.lt.1) iy=1 if (iy.gt.ghgyr2-ghgyr1+1) iy=ghgyr2-ghgyr1+1 xCH4=ghgam(3,iy)*H2ObyCH4 c If (AM_I_ROOT()) c & write(6,*) 'add in stratosphere: H2O gen. by CH4(ppm)=',xCH4 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) q(i,j,l)=q(i,j,l)+xCH4*dH2O(j,l,jmon)*byAM(l,i,j) #ifdef TRACERS_WATER C**** Add water to relevant tracers as well do n=1,ntm if (itime_tr0(n).le.itime) then select case (tr_wd_type(n)) case (nWater) ! water: add CH4-sourced water to tracers trm(i,j,l,n) = trm(i,j,l,n) + + tr_H2ObyCH4(n)*xCH4*dH2O(j,l,jmon)*dxyp(j) end select end if end do #endif aj(j,j_h2och4,:)=aj(j,j_h2och4,:)+ + xCH4*dH2O(j,l,jmon)*ftype(:,i,j) end do end do If (HAVE_NORTH_POLE) q(2:im,jm,l)=q(1,jm,l) If (HAVE_SOUTH_POLE) q(2:im, 1,l)=q(1, 1,l) #ifdef TRACERS_WATER do n=1,ntm If (HAVE_SOUTH_POLE) trm(2:im, 1,l,n)=trm(1, 1,l,n) If (HAVE_NORTH_POLE) trm(2:im,jm,l,n)=trm(1,jm,l,n) end do #endif end do end if RETURN END SUBROUTINE DAILY SUBROUTINE CHECKT (SUBR) 17,20 !@sum CHECKT Checks arrays for NaN/INF and reasonablness !@auth Original Development Team !@ver 1.0 C**** CHECKT IS TURNED ON BY SETTING QCHECK=.TRUE. IN NAMELIST C**** REMEMBER TO SET QCHECK BACK TO .FALSE. AFTER THE ERRORS ARE C**** CORRECTED. USE CONSTANT, only : tf USE MODEL_COM USE DYNAMICS, only : pk USE DOMAIN_DECOMP, only : grid, GET, AM_I_ROOT USE soil_drv, only : checke IMPLICIT NONE INTEGER I,J,L !@var SUBR identifies where CHECK was called from CHARACTER*6, INTENT(IN) :: SUBR c**** Extract domain decomposition info INTEGER :: J_0, J_1 CALL GET(grid, J_STRT = J_0, J_STOP = J_1) IF (QCHECK) THEN C**** Check all prog. arrays for Non-numbers CALL CHECK3(U,IM,JM,LM,SUBR,'u ') CALL CHECK3(V,IM,JM,LM,SUBR,'v ') CALL CHECK3(T,IM,JM,LM,SUBR,'t ') CALL CHECK3(Q,IM,JM,LM,SUBR,'q ') CALL CHECK3(P,IM,JM,1,SUBR,'p ') CALL CHECK3(WM,IM,JM,LM,SUBR,'wm ') DO J=J_0,J_1 DO I=1,IM IF (Q(I,J,1).gt.1d-1)print*,SUBR," Q BIG ",i,j,Q(I,J,1:LS1) IF (T(I,J,1)*PK(1,I,J)-TF.gt.50.) print*,SUBR," T BIG ",i,j * ,T(I,J,1:LS1)*PK(1:LS1,I,J)-TF END DO END DO DO L=1,LM DO J=J_0,J_1 DO I=1,IM IF (Q(I,J,L).lt.0.) then print*,"After ",SUBR," Q < 0 ",i,j,Q(I,J,L) call stop_model('Q<0 in CHECKT',255) END IF IF (WM(I,J,L).lt.0.) then print*,"After ",SUBR," WM < 0 ",i,j,WM(I,J,L) call stop_model('WM<0 in CHECKT',255) END IF END DO END DO END DO C**** Check PBL arrays CALL CHECKPBL(SUBR) C**** Check Ocean arrays CALL CHECKO(SUBR) C**** Check Ice arrays CALL CHECKI(SUBR) C**** Check Lake arrays CALL CHECKL(SUBR) C**** Check Earth arrays CALL CHECKE(SUBR) #if (defined TRACERS_ON) || (defined TRACERS_OCEAN) C**** check tracers CALL CHECKTR(SUBR) #endif END IF RETURN END SUBROUTINE CHECKT subroutine read_options( qcrestart, ifile ) 1,3 !@sum reads options from the command line (for now only one option) !@auth I. Aleinov !@ver 1.0 implicit none !@var qcrestart true if "-r" is present logical, intent(inout) :: qcrestart character(*),intent(inout) :: ifile character*80 arg,arg1 do call nextarg( arg, 1 ) if ( arg == "" ) exit ! end of args select case (arg) case ("-r") qcrestart = .true. case ("-i") call nextarg( arg1, 0 ) ifile=arg1 ! new options can be included here case default print *,'Unknown option specified: ', arg print *,'Aborting...' call stop_model("Unknown option on a command line",255) end select enddo return end subroutine read_options subroutine print_restart_info 1,11 !@sum prints timing information needed to restart the model !@auth I. Aleinov !@ver 1.0 USE MODEL_COM USE FILEMANAGER, only : openunit,closeunit use DOMAIN_DECOMP, only: AM_I_ROOT implicit none integer :: ItimeMax=-1, Itime1, Itime2, itm, ioerr1=-1, ioerr2=-1 integer :: iu_rsf call openunit(rsf_file_name(1),iu_rsf,.true.) call io_label(iu_rsf,Itime1,itm,ioread,ioerr1) call closeunit(iu_rsf) call openunit(rsf_file_name(2),iu_rsf,.true.) call io_label(iu_rsf,Itime2,itm,ioread,ioerr2) call closeunit(iu_rsf) if ( ioerr1==-1 ) ItimeMax = Itime1 if ( ioerr2==-1 ) ItimeMax = max( ItimeMax, Itime2 ) if ( Itime < 0 ) $ call stop_model("Could not read fort.1, fort.2",255) call getdte( & ItimeMax,Nday,Iyear1,Jyear,Jmon,Jday,Jdate,Jhour,amon) if (AM_I_ROOT()) & write(6,"('QCRESTART_DATA: ',I10,1X,I2,'-',I2.2,'-',I4.4)") & ItimeMax*24/Nday, Jmon, Jdate, Jyear return end subroutine print_restart_info