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