C $Header: /u/gcmpack/MITgcm/pkg/atm2d/forward_step_atm2d.F,v 1.27 2013/05/02 20:37:52 jmc Exp $ C $Name: $ #include "ctrparam.h" #ifdef OCEAN_3D # include "ATM2D_OPTIONS.h" #endif C SUBROUTINE FORWARD_STEP_ATM2D(iloop, myTime, myIter, myThid) C |==========================================================| C | Does time loop for one coupled period. The main loop | C | this is the MITGCM main loop OR a separate driver for | C | IGSM 2.2 | C \==========================================================/ IMPLICIT NONE #include "ATMSIZE.h" #include "DRIVER.h" #ifdef OCEAN_3D # include "SIZE.h" # include "EEPARAMS.h" # include "PARAMS.h" # include "ATM2D_VARS.h" #endif #ifdef NCEPWIND COMMON /SEED/JSEED,IFRST,NEXTN INTEGER JSEED,IFRST,NEXTN #endif C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == C iloop - loop counter for coupled period time steps (main time loop) C myIter - iteration counter for this thread (ocean times steps +nIter0) C myTime - time counter for this thread (ocean time, from starttTime) C myThid - thread number for this instance of the routine. INTEGER iloop REAL*8 myTime INTEGER myIter INTEGER myThid C === Local variables === INTEGER idyear ! year # of simulation, starting at year 1 INTEGER iyr ! year # of simulation, starting from specified inyear INTEGER inyr ! hours into the current year, end of coupled period INTEGER monid ! current month of the year INTEGER inday ! hour of the day, end of the coupled period INTEGER dayid ! day of the current month INTEGER j,mn,na,no !loop counters INTEGER jdofmhr(0:12) DATA jdofmhr/0,744,1416,2160,2880,3624,4344,5088, & 5832,6552,7296,8016,8760/ C i.e. 0,31*24,59*24,90*24,120*24,151*24,181*24, C 212*24,243*24,273*24,304*24,334*24,365*24 #ifdef CPL_TEM INTEGER ndmonth(12) DATA ndmonth/31,28,31,30,31,30,31,31,30,31,30,31/ REAL*4 totup, aduptt REAL*8 tcumn #endif #if (defined CPL_TEM) (defined CPL_OCEANCO2) REAL*8 nepan REAL*8 ocuptp #endif #ifdef OCEAN_3D INTEGER iloop_ocn, i _RL qPrcRn(1-OLx:sNx+OLx,1-OLy:sNy+OLy) # ifdef NCEPWIND REAL*8 RAND CHARACTER *4 ncep_yr # endif #endif #ifdef DATA4TEM CHARACTER *40 f4tem,f4clm,f24tem,f34tem CHARACTER *4 cfile CHARACTER *8 f14tem,f14clm character *9 f124tem,f134tem f14tem='data4tem' f14clm='data4clm' f124tem='data24tem' f134tem='data34tem' nfile=1 #endif C print *,'*** Top of forwrdstep_atm',iloop,myTime,myIter idyear= int((iloop-1)*dtcouple/365.0/24.0) + 1 iyr= idyear + startYear -1 inyr = mod(iloop*dtcouple, 365*24) IF (inyr .EQ. 0) inyr=jdofmhr(12) DO mn=1,12 IF ((inyr.GT.jdofmhr(mn-1)).AND.(inyr.LE.jdofmhr(mn))) monid=mn ENDDO inday= mod(iloop*dtcouple, 24) dayid= int((inyr-dtcouple-jdofmhr(monid-1))/24.0) +1 C print *,'*** idyear,iyr,inyr,monid,inday,dayid', C & idyear,iyr,inyr,monid,inday,dayid IF (inyr.EQ.dtcouple) THEN !do this block at start of new year PRINT *,'*** Starting a new year' #ifdef NCEPWIND WRITE(ncep_yr,'(I4)') (NINT(RAND()*60.0+0.5) + 1947) PRINT *,'Using NCEP wind variations from year: ',ncep_yr OPEN(6007, & FILE='ncep_taux_variations_'//ncep_yr//'.bin',STATUS='old', & ACCESS='direct', RECL=4*sNx*sNy, & FORM='unformatted') OPEN(6008, & FILE='ncep_tauy_variations_'//ncep_yr//'.bin',STATUS='old', & ACCESS='direct', RECL=4*sNx*sNy, & FORM='unformatted') OPEN(6009, & FILE='ncep_speed_variations_'//ncep_yr//'.bin',STATUS='old', & ACCESS='direct', RECL=4*sNx*sNy, & FORM='unformatted') ncep_counter = 1 #endif #ifdef DATA4TEM IF (nfile.gt.1)THEN CLOSE(935) CLOSE(937) CLOSE(938) CLOSE(939) ENDIF IF(iyr.gt.1000) THEN nfile=iyr ELSE nfile=1000+iyr ENDIF WRITE (cfile,'i4'),nfile f4tem=f14tem//cfile f4clm=f14clm//cfile f24tem=f124tem//cfile f34tem=f134tem//cfile OPEN(935,file=f4clm,form='unformatted',status='new') OPEN(937,file=f4tem,form='unformatted',status='new') OPEN(938,file=f24tem,form='unformatted',status='new') OPEN(939,file=f34tem,form='unformatted',status='new') nfile=nfile+1 #endif #ifdef CPL_TEM nepan=0.0 ch4ann=0.0 n2oann=0.0 xco2ann=0.0 #endif #ifdef CPL_OCEANCO2 temuptann=0. DO j=1,jm0 co24ocnan(j)=0.0 ENDDO # ifdef ML_2D call KVCARBON(iyr) # endif #endif #ifdef CPL_TEM DO j=1,jm0 antemnep(j)=0. ENDDO # ifndef CPL_CHEM CALL ROBSO3(iyr) # endif C For land use CALL UPDATELCLUC(idyear) #endif #ifdef CPL_CHEM print *,' Before eppaemission' CALL EPPAEMISSION (iyr) #endif ENDIF !end block done at year-start IF (inyr.EQ.jdofmhr(monid-1)+dtcouple) THEN !do this block month start PRINT *,'***Starting a new month' #ifdef CPL_TEM CALL ZCLIMATE2TEM #endif #ifdef CPL_OCEANCO2 ocumn=0.0 DO j=1,jm0 fluxco2mn(j)=0.0 ENDDO #endif #ifdef OCEAN_3D new_mon= .TRUE. #endif ENDIF !end block at start of the month C C------------------- Top of Coupled Period Loop -------------------------- C #ifdef OCEAN_3D # ifdef NCEPWIND C Read in the next timestep of ncep wind stress variations C PRINT *,'*** Read in next NCEP record' READ(6007, REC=ncep_counter), fu_ncep READ(6008, REC=ncep_counter), fv_ncep READ(6009, REC=ncep_counter), fs_ncep ncep_counter=ncep_counter+1 # endif # ifdef ATM2D_MPI_ON CALL CPL_RECV_OCN_FIELDS # endif CALL GET_OCNVARS( myTime, myIter, myThid) IF ( (iloop.NE.1).OR. (iloop.EQ.1.AND. & (startTime.NE.baseTime .OR. nIter0.NE.0)) ) THEN C don't run the ice growth/melt on step 1 if "cold" start DO j = 1-OLy, sNy+OLy DO i = 1-OLx, sNx+OLx qPrcRn(i,j) = 0. ENDDO ENDDO CALL THSICE_STEP_FWD( 1, 1, 1, sNx, 1, sNy, I pass_prcAtm, snowPrc, qPrcRn, & myTime, myIter, myThid ) CALL THSICE_AVE( 1,1, myTime, myIter, myThid ) ENDIF CALL CALC_ZONAL_MEANS(.TRUE.,myThid) CALL PUT_OCNVARS(myTime,myIter,myThid) CALL SUM_YR_END_DIAGS(myTime,myIter,myThid) # ifdef ATM2D_MPI_ON CALL CPL_SEND_OCN_FIELDS # endif #endif C PRINT *,'Top of ncall_atm Loop' DO na=1,ncall_atm !loop for atmos forward time steps CALL ATMOSPHERE(dtatm,monid) #ifdef OCEAN_3D CALL ATM2OCN_MAIN(iloop, na, monid, myIter, myThid) CALL SUM_OCN_FLUXES(myThid) CALL PASS_THSICE_FLUXES(myThid) CALL THSICE_IMPL_TEMP(netSW, sFlx, dTsurf, 1,1, & myTime, myIter, myThid) CALL SUM_THSICE_OUT(myThid) CALL CALC_ZONAL_MEANS(.FALSE.,myThid) !just mean Tsrf recalculated #endif ENDDO ! ncall_atm loop C PRINT *,'Top of ncall_ocean Loop' DO no=1,ncall_ocean !loop for each ocean forward step #ifdef OCEAN_3D iloop_ocn = nint((iloop-1)*dtcouplo/deltaTClock) + no # ifndef ATM2D_MPI_ON CALL FORWARD_STEP(iloop_ocn, myTime, myIter, myThid ) else myIter = nIter0 + iloop_ocn myTime = startTime + deltaTClock *float (iloop_ocn) CALL DO_THE_MODEL_IO( .FALSE., myTime, myIter, myThid ) CALL DO_WRITE_PICKUP( .FALSE., myTime, myIter, myThid ) # endif #endif #ifdef ML_2D CALL OCEAN_ML(dtocn*3600.,dtatm*3600.) #endif ENDDO ! ncall_ocean loop C Start of code done at the end of every coupled period #ifdef OCEAN_3D CALL NORM_OCN_FLUXES(myThid) CALL ATM2D_WRITE_PICKUP(.FALSE., myTime, myIter, myThid) #endif C C--------------------- End of coupled period loop -------------------- C IF (inday.EQ.0) THEN !do this block if end-of-day C PRINT *,'***block at end of day' ENDIF !end block end-of-day IF (inyr.EQ.jdofmhr(monid)) THEN !do block if month-end PRINT *,'***end of month reached' #ifdef CLM # ifdef CPL_TEM CALL CLIMATE2TEM(monid,ndmonth(monid)) c PRINT *,'From driver before call tem',' idyear=',idyear CALL TEM(idyear,monid-1) CALL TEM2CLIMATE(idyear,monid-1) ch4mn=0.0 n2omn=0.0 nepmn=0.0 DO j=1,jm0 ch4mn=ch4mn+temch4(j) n2omn=n2omn+temn2o(j) nepmn=nepmn+temco2(j) ENDDO # ifdef CPL_NEM PRINT *,'Month=',monid PRINT *,'CH4=',ch4mn/1.e9,' N2O=',n2omn/1.e9 OPEN(277,ACCESS='APPEND',FILE=fnememiss,form='unformatted' & ,STATUS='old') WRITE (277) iyr,monid,ch4mn,n2omn,nepmn, & temch4,temn2o,temco2 CLOSE(277) # endif DO j=1,jm0 temnep(monid,j)=temco2(j) ENDDO DO j=1,jm0 antemnep(j)=antemnep(j)+temnep(monid,j) nepan=nepan+temnep(monid,j) ch4ann=ch4ann+temch4(j) n2oann=n2oann+temn2o(j) ENDDO # endif #endif #ifdef OCEAN_3D CALL MONTH_END_DIAGS( monid, myTime, myIter, myThid) #endif #ifdef CPL_OCEANCO2 IF (monid.EQ.12) THEN ocupt=ocupt*12.e-15 C 12.e-15 from moles to Gt carbon ocuptp=ocupt ocupt=0.0 ENDIF #endif #ifdef IPCC_EMI PRINT *,'Month=',monid nepmn=nepmn/1000. C nepmn is converted from Tg/mn to Gt/mn of C ocumn=ocumn*12.e-15 C ocumn is converted from mole(C) to Gt (C) C tnow= jyear + (jday-.5)/365. C CALL emissipcc(tnow,nepmn,ocumn,CO2,xco2ann) print *,nepmn,ocumn,xco2ann C ch4mn is in kg/mn of CH4 C nepmn is in Gt/mn of C tcumn=nepmn-ch4mn*12./16.*1.e-12 print *,'tcumn,ocumn,xco2ann' print *,tcumn,ocumn,xco2ann CALL EMISSIPCC_MN(tcumn,ocumn,xco2ann) C CALL emissipcc_mn(nepmn,ocumn,xco2ann) #endif ENDIF !end block done at month-end IF (inyr .EQ. jdofmhr(12)) THEN ! do this block at year-end PRINT *,'***end of year reached' #ifdef CPL_TEM nepan=nepan/1000. IF (iyr.ge.1981.and.iyr.le.1990) THEN PRINT *,'Uptake avegaging year=',iyr nepav=nepav+nepan aocuav=aocuav+OCUPTP IF (iyr.eq.1990) THEN nepav=nepav/10. aocuav=aocuav/10. totup=nepav+aocuav aduptt=4.1-totup PRINT *,' Carbon uptake for spinup' PRINT *,' totup=',totup,' aocuav=',aocuav PRINT *,' nepav=',nepav,' aduptt=',aduptt ENDIF ENDIF IF (iyr.eq.endYear) THEN C NEM emissions and NEP for start of climate-chemistry run adupt=aduptt CALL WR_RSTRT_NEM ENDIF #endif #ifdef ML_2D C Data for the restart of the 2D ML model CALL WRRSTRT_OCEAN #endif #ifdef OCEAN_3D IF ((mod(iyr,taveDump).EQ.0).AND.(idyear.GE.taveDump)) THEN CALL TAVE_END_DIAGS( taveDump, myTime, myIter, myThid) ELSEIF (mod(iyr,taveDump).EQ.0) THEN CALL TAVE_END_DIAGS( idyear, myTime, myIter, myThid) ENDIF CALL YR_END_DIAGS(iyr,myTime,myIter,myThid) C If necessary, next line can be moved outside OCEAN_3D for IGSM2.2 cleanups IF (iloop.EQ.nTimeSteps) CALL ATM2D_FINISH(myThid) # ifdef NCEPWIND OPEN(unit=334,file='rand_state_new.dat',status='replace') WRITE(334,*) JSEED,IFRST,NEXTN CLOSE(334) # endif #endif #if (defined CPL_TEM) (defined CPL_OCEANCO2) PRINT '(a6,i6,2(a5,f10.4))','Year=',iyr, & ' NEP=',nepan,' OCU=',OCUPTP # ifdef IPCC_EMI PRINT '(a6,i6,(a5,f10.4))','Year=',iyr, & ' CO2AN=',xco2ann/12. CALL EMISSIPCC_YR # endif # ifdef CPL_NEM PRINT *,' CH4=',ch4ann,' N2O=',n2oann # endif OPEN(333,ACCESS='APPEND',FILE=caruptfile,STATUS='old') # ifndef CPL_TEM C For ocean carbon model only WRITE(333,'(i7,3e15.5)')iyr,ocuptp else # ifndef CPL_OCEANCO2 C For ocean TEM only WRITE(333,'(i7,3e15.5)')iyr,nepan,nepan-1.e-12*ch4ann*12./16. else C For ocean both TEM OCM WRITE(333,'(i7,3e15.5)')iyr,nepan,nepan-1.e-12*ch4ann*12./16. & ,ocuptp # endif # endif CLOSE(333) # if (defined CPL_OCEANCO2) (defined ML_2D) WRITE(602)iyr CALL WRGARY CALL ZEROGARY # endif #endif #ifdef CPL_OCEANCO2 DO j=1,jm0 # ifdef OCEAN_3D co24ocnan(j)=co24ocnan(j)*dtatm/24.0/365.0 else co24ocnan(j)=co24ocnan(j)/365.0 # endif ENDDO PRINT *,' CO2 for ocean model',' ncallatm=',ncall_atm PRINT '(12f7.1,/,2(11f7.1,/),12f7.1)',co24ocnan #endif #ifdef CPL_CHEM PRINT *,' TEMUPTANN=',temuptann,' TOTAL UPTAKE=' & ,ocuptp+temuptann #endif ENDIF !year-end block RETURN END
#ifdef NCEPWIND REAL*8 FUNCTION RAND() C C This function returns a pseudo-random number for each invocation. C It is a FORTRAN 77 adaptation of the "Integer Version 2" minimal C standard number generator whose Pascal code appears in the article: C C Park, Steven K. and Miller, Keith W., "Random Number Generators: C Good Ones are Hard to Find", Communications of the ACM, C October, 1988. C PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773, + MOMDMP=2836) C COMMON /SEED/JSEED,IFRST,NEXTN INTEGER JSEED,IFRST,NEXTN INTEGER HVLUE, LVLUE, TESTV C IF (IFRST .EQ. 0) THEN NEXTN = JSEED IFRST = 1 ENDIF C HVLUE = NEXTN / MOBYMP LVLUE = MOD(NEXTN, MOBYMP) TESTV = MPLIER*LVLUE - MOMDMP*HVLUE IF (TESTV .GT. 0) THEN NEXTN = TESTV ELSE NEXTN = TESTV + MODLUS ENDIF RAND = REAL(NEXTN)/REAL(MODLUS) C RETURN END
#endif