1 |
C $Header: /u/gcmpack/MITgcm/pkg/atm2d/forward_step_atm2d.F,v 1.26 2013/04/04 17:06:37 jmc Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "ctrparam.h" |
5 |
#ifdef OCEAN_3D |
6 |
# include "ATM2D_OPTIONS.h" |
7 |
#endif |
8 |
C |
9 |
SUBROUTINE FORWARD_STEP_ATM2D(iloop, myTime, myIter, myThid) |
10 |
C |==========================================================| |
11 |
C | Does time loop for one coupled period. The main loop | |
12 |
C | this is the MITGCM main loop OR a separate driver for | |
13 |
C | IGSM 2.2 | |
14 |
C \==========================================================/ |
15 |
IMPLICIT NONE |
16 |
|
17 |
#include "ATMSIZE.h" |
18 |
#include "DRIVER.h" |
19 |
|
20 |
#ifdef OCEAN_3D |
21 |
# include "SIZE.h" |
22 |
# include "EEPARAMS.h" |
23 |
# include "PARAMS.h" |
24 |
# include "ATM2D_VARS.h" |
25 |
#endif |
26 |
|
27 |
#ifdef NCEPWIND |
28 |
COMMON /SEED/JSEED,IFRST,NEXTN |
29 |
INTEGER JSEED,IFRST,NEXTN |
30 |
#endif |
31 |
|
32 |
C !INPUT/OUTPUT PARAMETERS: |
33 |
C == Routine arguments == |
34 |
C iloop - loop counter for coupled period time steps (main time loop) |
35 |
C myIter - iteration counter for this thread (ocean times steps +nIter0) |
36 |
C myTime - time counter for this thread (ocean time, from starttTime) |
37 |
C myThid - thread number for this instance of the routine. |
38 |
INTEGER iloop |
39 |
REAL*8 myTime |
40 |
INTEGER myIter |
41 |
INTEGER myThid |
42 |
|
43 |
C === Local variables === |
44 |
INTEGER idyear ! year # of simulation, starting at year 1 |
45 |
INTEGER iyr ! year # of simulation, starting from specified inyear |
46 |
INTEGER inyr ! hours into the current year, end of coupled period |
47 |
INTEGER monid ! current month of the year |
48 |
INTEGER inday ! hour of the day, end of the coupled period |
49 |
INTEGER dayid ! day of the current month |
50 |
INTEGER j,mn,na,no !loop counters |
51 |
INTEGER jdofmhr(0:12) |
52 |
DATA jdofmhr/0,744,1416,2160,2880,3624,4344,5088, |
53 |
& 5832,6552,7296,8016,8760/ |
54 |
C i.e. 0,31*24,59*24,90*24,120*24,151*24,181*24, |
55 |
C 212*24,243*24,273*24,304*24,334*24,365*24 |
56 |
#ifdef CPL_TEM |
57 |
INTEGER ndmonth(12) |
58 |
DATA ndmonth/31,28,31,30,31,30,31,31,30,31,30,31/ |
59 |
REAL*4 totup, aduptt |
60 |
REAL*8 tcumn |
61 |
#endif |
62 |
#if (defined CPL_TEM) || (defined CPL_OCEANCO2) |
63 |
REAL*8 nepan |
64 |
REAL*8 ocuptp |
65 |
#endif |
66 |
#ifdef OCEAN_3D |
67 |
INTEGER iloop_ocn, i |
68 |
_RL qPrcRn(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
69 |
# ifdef NCEPWIND |
70 |
REAL*8 RAND |
71 |
CHARACTER *4 ncep_yr |
72 |
# endif |
73 |
#endif |
74 |
#ifdef DATA4TEM |
75 |
CHARACTER *40 f4tem,f4clm,f24tem,f34tem |
76 |
CHARACTER *4 cfile |
77 |
CHARACTER *8 f14tem,f14clm |
78 |
character *9 f124tem,f134tem |
79 |
f14tem='data4tem' |
80 |
f14clm='data4clm' |
81 |
f124tem='data24tem' |
82 |
f134tem='data34tem' |
83 |
nfile=1 |
84 |
#endif |
85 |
|
86 |
C print *,'*** Top of forwrdstep_atm',iloop,myTime,myIter |
87 |
idyear= int((iloop-1)*dtcouple/365.0/24.0) + 1 |
88 |
iyr= idyear + startYear -1 |
89 |
inyr = mod(iloop*dtcouple, 365*24) |
90 |
IF (inyr .EQ. 0) inyr=jdofmhr(12) |
91 |
DO mn=1,12 |
92 |
IF ((inyr.GT.jdofmhr(mn-1)).AND.(inyr.LE.jdofmhr(mn))) monid=mn |
93 |
ENDDO |
94 |
inday= mod(iloop*dtcouple, 24) |
95 |
dayid= int((inyr-dtcouple-jdofmhr(monid-1))/24.0) +1 |
96 |
C print *,'*** idyear,iyr,inyr,monid,inday,dayid', |
97 |
C & idyear,iyr,inyr,monid,inday,dayid |
98 |
|
99 |
IF (inyr.EQ.dtcouple) THEN !do this block at start of new year |
100 |
PRINT *,'*** Starting a new year' |
101 |
#ifdef NCEPWIND |
102 |
WRITE(ncep_yr,'(I4)') (NINT(RAND()*60.0+0.5) + 1947) |
103 |
PRINT *,'Using NCEP wind variations from year: ',ncep_yr |
104 |
OPEN(6007, |
105 |
& FILE='ncep_taux_variations_'//ncep_yr//'.bin',STATUS='old', |
106 |
& ACCESS='direct', RECL=4*sNx*sNy, |
107 |
& FORM='unformatted') |
108 |
OPEN(6008, |
109 |
& FILE='ncep_tauy_variations_'//ncep_yr//'.bin',STATUS='old', |
110 |
& ACCESS='direct', RECL=4*sNx*sNy, |
111 |
& FORM='unformatted') |
112 |
OPEN(6009, |
113 |
& FILE='ncep_speed_variations_'//ncep_yr//'.bin',STATUS='old', |
114 |
& ACCESS='direct', RECL=4*sNx*sNy, |
115 |
& FORM='unformatted') |
116 |
ncep_counter = 1 |
117 |
#endif |
118 |
#ifdef DATA4TEM |
119 |
IF (nfile.gt.1)THEN |
120 |
CLOSE(935) |
121 |
CLOSE(937) |
122 |
CLOSE(938) |
123 |
CLOSE(939) |
124 |
ENDIF |
125 |
IF(iyr.gt.1000) THEN |
126 |
nfile=iyr |
127 |
ELSE |
128 |
nfile=1000+iyr |
129 |
ENDIF |
130 |
WRITE (cfile,'i4'),nfile |
131 |
f4tem=f14tem//cfile |
132 |
f4clm=f14clm//cfile |
133 |
f24tem=f124tem//cfile |
134 |
f34tem=f134tem//cfile |
135 |
OPEN(935,file=f4clm,form='unformatted',status='new') |
136 |
OPEN(937,file=f4tem,form='unformatted',status='new') |
137 |
OPEN(938,file=f24tem,form='unformatted',status='new') |
138 |
OPEN(939,file=f34tem,form='unformatted',status='new') |
139 |
nfile=nfile+1 |
140 |
#endif |
141 |
#ifdef CPL_TEM |
142 |
nepan=0.0 |
143 |
ch4ann=0.0 |
144 |
n2oann=0.0 |
145 |
xco2ann=0.0 |
146 |
#endif |
147 |
#ifdef CPL_OCEANCO2 |
148 |
temuptann=0. |
149 |
DO j=1,jm0 |
150 |
co24ocnan(j)=0.0 |
151 |
ENDDO |
152 |
# ifdef ML_2D |
153 |
call kvcarbon(iyr) |
154 |
# endif |
155 |
#endif |
156 |
#ifdef CPL_TEM |
157 |
DO j=1,jm0 |
158 |
antemnep(j)=0. |
159 |
ENDDO |
160 |
# ifndef CPL_CHEM |
161 |
CALL robso3(iyr) |
162 |
# endif |
163 |
C For land use |
164 |
CALL updatelcluc(idyear) |
165 |
#endif |
166 |
#ifdef CPL_CHEM |
167 |
print *,' Before eppaemission' |
168 |
CALL eppaemission (iyr) |
169 |
#endif |
170 |
ENDIF !end block done at year-start |
171 |
|
172 |
IF (inyr.EQ.jdofmhr(monid-1)+dtcouple) THEN !do this block month start |
173 |
PRINT *,'***Starting a new month' |
174 |
#ifdef CPL_TEM |
175 |
CALL zclimate2tem |
176 |
#endif |
177 |
#ifdef CPL_OCEANCO2 |
178 |
ocumn=0.0 |
179 |
DO j=1,jm0 |
180 |
fluxco2mn(j)=0.0 |
181 |
ENDDO |
182 |
#endif |
183 |
#ifdef OCEAN_3D |
184 |
new_mon= .TRUE. |
185 |
#endif |
186 |
ENDIF !end block at start of the month |
187 |
C |
188 |
C------------------- Top of Coupled Period Loop -------------------------- |
189 |
C |
190 |
|
191 |
#ifdef OCEAN_3D |
192 |
# ifdef NCEPWIND |
193 |
C Read in the next timestep of ncep wind stress variations |
194 |
C PRINT *,'*** Read in next NCEP record' |
195 |
READ(6007, REC=ncep_counter), fu_ncep |
196 |
READ(6008, REC=ncep_counter), fv_ncep |
197 |
READ(6009, REC=ncep_counter), fs_ncep |
198 |
ncep_counter=ncep_counter+1 |
199 |
# endif |
200 |
# ifdef ATM2D_MPI_ON |
201 |
CALL CPL_RECV_OCN_FIELDS |
202 |
# endif |
203 |
CALL GET_OCNVARS( myTime, myIter, myThid) |
204 |
IF ( (iloop.NE.1).OR. (iloop.EQ.1.AND. |
205 |
& (startTime.NE.baseTime .OR. nIter0.NE.0)) ) THEN |
206 |
C don't run the ice growth/melt on step 1 if "cold" start |
207 |
DO j = 1-OLy, sNy+OLy |
208 |
DO i = 1-OLx, sNx+OLx |
209 |
qPrcRn(i,j) = 0. |
210 |
ENDDO |
211 |
ENDDO |
212 |
CALL THSICE_STEP_FWD( 1, 1, 1, sNx, 1, sNy, |
213 |
I pass_prcAtm, snowPrc, qPrcRn, |
214 |
& myTime, myIter, myThid ) |
215 |
CALL THSICE_AVE( 1,1, myTime, myIter, myThid ) |
216 |
ENDIF |
217 |
CALL CALC_ZONAL_MEANS(.TRUE.,myThid) |
218 |
CALL PUT_OCNVARS(myTime,myIter,myThid) |
219 |
CALL SUM_YR_END_DIAGS(myTime,myIter,myThid) |
220 |
# ifdef ATM2D_MPI_ON |
221 |
CALL CPL_SEND_OCN_FIELDS |
222 |
# endif |
223 |
#endif |
224 |
|
225 |
C PRINT *,'Top of ncall_atm Loop' |
226 |
DO na=1,ncall_atm !loop for atmos forward time steps |
227 |
CALL atmosphere(dtatm,monid) |
228 |
#ifdef OCEAN_3D |
229 |
CALL ATM2OCN_MAIN(iloop, na, monid, myIter, myThid) |
230 |
CALL SUM_OCN_FLUXES(myThid) |
231 |
CALL PASS_THSICE_FLUXES(myThid) |
232 |
CALL THSICE_IMPL_TEMP(netSW, sFlx, dTsurf, 1,1, |
233 |
& myTime, myIter, myThid) |
234 |
CALL SUM_THSICE_OUT(myThid) |
235 |
CALL CALC_ZONAL_MEANS(.FALSE.,myThid) !just mean Tsrf recalculated |
236 |
#endif |
237 |
ENDDO ! ncall_atm loop |
238 |
|
239 |
C PRINT *,'Top of ncall_ocean Loop' |
240 |
DO no=1,ncall_ocean !loop for each ocean forward step |
241 |
|
242 |
#ifdef OCEAN_3D |
243 |
iloop_ocn = nint((iloop-1)*dtcouplo/deltaTClock) + no |
244 |
# ifndef ATM2D_MPI_ON |
245 |
CALL FORWARD_STEP(iloop_ocn, myTime, myIter, myThid ) |
246 |
# else |
247 |
myIter = nIter0 + iloop_ocn |
248 |
myTime = startTime + deltaTClock *float (iloop_ocn) |
249 |
CALL DO_THE_MODEL_IO( .FALSE., myTime, myIter, myThid ) |
250 |
CALL DO_WRITE_PICKUP( .FALSE., myTime, myIter, myThid ) |
251 |
# endif |
252 |
#endif |
253 |
#ifdef ML_2D |
254 |
CALL ocean_ml(dtocn*3600.,dtatm*3600.) |
255 |
#endif |
256 |
|
257 |
ENDDO ! ncall_ocean loop |
258 |
|
259 |
C Start of code done at the end of every coupled period |
260 |
|
261 |
#ifdef OCEAN_3D |
262 |
CALL NORM_OCN_FLUXES(myThid) |
263 |
CALL ATM2D_WRITE_PICKUP(.FALSE., myTime, myIter, myThid) |
264 |
#endif |
265 |
|
266 |
C |
267 |
C--------------------- End of coupled period loop -------------------- |
268 |
C |
269 |
IF (inday.EQ.0) THEN !do this block if end-of-day |
270 |
C PRINT *,'***block at end of day' |
271 |
ENDIF !end block end-of-day |
272 |
|
273 |
IF (inyr.EQ.jdofmhr(monid)) THEN !do block if month-end |
274 |
PRINT *,'***end of month reached' |
275 |
#ifdef CLM |
276 |
# ifdef CPL_TEM |
277 |
CALL climate2tem(monid,ndmonth(monid)) |
278 |
c PRINT *,'From driver before call tem',' idyear=',idyear |
279 |
CALL tem(idyear,monid-1) |
280 |
CALL tem2climate(idyear,monid-1) |
281 |
ch4mn=0.0 |
282 |
n2omn=0.0 |
283 |
nepmn=0.0 |
284 |
DO j=1,jm0 |
285 |
ch4mn=ch4mn+temch4(j) |
286 |
n2omn=n2omn+temn2o(j) |
287 |
nepmn=nepmn+temco2(j) |
288 |
ENDDO |
289 |
# ifdef CPL_NEM |
290 |
PRINT *,'Month=',monid |
291 |
PRINT *,'CH4=',ch4mn/1.e9,' N2O=',n2omn/1.e9 |
292 |
OPEN(277,ACCESS='APPEND',FILE=fnememiss,form='unformatted' |
293 |
& ,STATUS='old') |
294 |
WRITE (277) iyr,monid,ch4mn,n2omn,nepmn, |
295 |
& temch4,temn2o,temco2 |
296 |
CLOSE(277) |
297 |
# endif |
298 |
DO j=1,jm0 |
299 |
temnep(monid,j)=temco2(j) |
300 |
ENDDO |
301 |
DO j=1,jm0 |
302 |
antemnep(j)=antemnep(j)+temnep(monid,j) |
303 |
nepan=nepan+temnep(monid,j) |
304 |
ch4ann=ch4ann+temch4(j) |
305 |
n2oann=n2oann+temn2o(j) |
306 |
ENDDO |
307 |
|
308 |
# endif |
309 |
#endif |
310 |
|
311 |
#ifdef OCEAN_3D |
312 |
CALL MONTH_END_DIAGS( monid, myTime, myIter, myThid) |
313 |
#endif |
314 |
|
315 |
#ifdef CPL_OCEANCO2 |
316 |
IF (monid.EQ.12) THEN |
317 |
ocupt=ocupt*12.e-15 |
318 |
C 12.e-15 from moles to Gt carbon |
319 |
ocuptp=ocupt |
320 |
ocupt=0.0 |
321 |
ENDIF |
322 |
#endif |
323 |
|
324 |
#ifdef IPCC_EMI |
325 |
PRINT *,'Month=',monid |
326 |
nepmn=nepmn/1000. |
327 |
C nepmn is converted from Tg/mn to Gt/mn of C |
328 |
ocumn=ocumn*12.e-15 |
329 |
C ocumn is converted from mole(C) to Gt (C) |
330 |
C tnow= jyear + (jday-.5)/365. |
331 |
C CALL emissipcc(tnow,nepmn,ocumn,CO2,xco2ann) |
332 |
print *,nepmn,ocumn,xco2ann |
333 |
C ch4mn is in kg/mn of CH4 |
334 |
C nepmn is in Gt/mn of C |
335 |
tcumn=nepmn-ch4mn*12./16.*1.e-12 |
336 |
print *,'tcumn,ocumn,xco2ann' |
337 |
print *,tcumn,ocumn,xco2ann |
338 |
CALL emissipcc_mn(tcumn,ocumn,xco2ann) |
339 |
C CALL emissipcc_mn(nepmn,ocumn,xco2ann) |
340 |
#endif |
341 |
ENDIF !end block done at month-end |
342 |
|
343 |
IF (inyr .EQ. jdofmhr(12)) THEN ! do this block at year-end |
344 |
PRINT *,'***end of year reached' |
345 |
#ifdef CPL_TEM |
346 |
nepan=nepan/1000. |
347 |
IF (iyr.ge.1981.and.iyr.le.1990) THEN |
348 |
PRINT *,'Uptake avegaging year=',iyr |
349 |
nepav=nepav+nepan |
350 |
aocuav=aocuav+OCUPTP |
351 |
IF (iyr.eq.1990) THEN |
352 |
nepav=nepav/10. |
353 |
aocuav=aocuav/10. |
354 |
totup=nepav+aocuav |
355 |
aduptt=4.1-totup |
356 |
PRINT *,' Carbon uptake for spinup' |
357 |
PRINT *,' totup=',totup,' aocuav=',aocuav |
358 |
PRINT *,' nepav=',nepav,' aduptt=',aduptt |
359 |
ENDIF |
360 |
ENDIF |
361 |
|
362 |
IF (iyr.eq.endYear) THEN |
363 |
C NEM emissions and NEP for start of climate-chemistry run |
364 |
adupt=aduptt |
365 |
CALL wr_rstrt_nem |
366 |
ENDIF |
367 |
|
368 |
#endif |
369 |
|
370 |
#ifdef ML_2D |
371 |
C Data for the restart of the 2D ML model |
372 |
CALL wrrstrt_ocean |
373 |
#endif |
374 |
|
375 |
#ifdef OCEAN_3D |
376 |
IF ((mod(iyr,taveDump).EQ.0).AND.(idyear.GE.taveDump)) THEN |
377 |
CALL TAVE_END_DIAGS( taveDump, myTime, myIter, myThid) |
378 |
ELSEIF (mod(iyr,taveDump).EQ.0) THEN |
379 |
CALL TAVE_END_DIAGS( idyear, myTime, myIter, myThid) |
380 |
ENDIF |
381 |
CALL YR_END_DIAGS(iyr,myTime,myIter,myThid) |
382 |
C If necessary, next line can be moved outside OCEAN_3D for IGSM2.2 cleanups |
383 |
IF (iloop.EQ.nTimeSteps) CALL ATM2D_FINISH(myThid) |
384 |
# ifdef NCEPWIND |
385 |
OPEN(unit=334,file='rand_state_new.dat',status='replace') |
386 |
WRITE(334,*) JSEED,IFRST,NEXTN |
387 |
CLOSE(334) |
388 |
# endif |
389 |
#endif |
390 |
|
391 |
#if (defined CPL_TEM) || (defined CPL_OCEANCO2) |
392 |
PRINT '(a6,i6,2(a5,f10.4))','Year=',iyr, |
393 |
& ' NEP=',nepan,' OCU=',OCUPTP |
394 |
# ifdef IPCC_EMI |
395 |
PRINT '(a6,i6,(a5,f10.4))','Year=',iyr, |
396 |
& ' CO2AN=',xco2ann/12. |
397 |
CALL emissipcc_yr |
398 |
# endif |
399 |
|
400 |
# ifdef CPL_NEM |
401 |
PRINT *,' CH4=',ch4ann,' N2O=',n2oann |
402 |
# endif |
403 |
OPEN(333,ACCESS='APPEND',FILE=caruptfile,STATUS='old') |
404 |
# ifndef CPL_TEM |
405 |
C For ocean carbon model only |
406 |
WRITE(333,'(i7,3e15.5)')iyr,ocuptp |
407 |
# else |
408 |
# ifndef CPL_OCEANCO2 |
409 |
C For ocean TEM only |
410 |
WRITE(333,'(i7,3e15.5)')iyr,nepan,nepan-1.e-12*ch4ann*12./16. |
411 |
# else |
412 |
C For ocean both TEM OCM |
413 |
WRITE(333,'(i7,3e15.5)')iyr,nepan,nepan-1.e-12*ch4ann*12./16. |
414 |
& ,ocuptp |
415 |
# endif |
416 |
# endif |
417 |
CLOSE(333) |
418 |
# if (defined CPL_OCEANCO2) && (defined ML_2D) |
419 |
WRITE(602)iyr |
420 |
CALL wrgary |
421 |
CALL zerogary |
422 |
# endif |
423 |
#endif |
424 |
|
425 |
#ifdef CPL_OCEANCO2 |
426 |
DO j=1,jm0 |
427 |
# ifdef OCEAN_3D |
428 |
co24ocnan(j)=co24ocnan(j)*dtatm/24.0/365.0 |
429 |
# else |
430 |
co24ocnan(j)=co24ocnan(j)/365.0 |
431 |
# endif |
432 |
ENDDO |
433 |
PRINT *,' CO2 for ocean model',' ncallatm=',ncall_atm |
434 |
PRINT '(12f7.1,/,2(11f7.1,/),12f7.1)',co24ocnan |
435 |
#endif |
436 |
|
437 |
#ifdef CPL_CHEM |
438 |
PRINT *,' TEMUPTANN=',temuptann,' TOTAL UPTAKE=' |
439 |
& ,ocuptp+temuptann |
440 |
#endif |
441 |
ENDIF !year-end block |
442 |
|
443 |
RETURN |
444 |
END |
445 |
|
446 |
#ifdef NCEPWIND |
447 |
REAL*8 FUNCTION RAND() |
448 |
C |
449 |
C This function returns a pseudo-random number for each invocation. |
450 |
C It is a FORTRAN 77 adaptation of the "Integer Version 2" minimal |
451 |
C standard number generator whose Pascal code appears in the article: |
452 |
C |
453 |
C Park, Steven K. and Miller, Keith W., "Random Number Generators: |
454 |
C Good Ones are Hard to Find", Communications of the ACM, |
455 |
C October, 1988. |
456 |
C |
457 |
PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773, |
458 |
+ MOMDMP=2836) |
459 |
C |
460 |
COMMON /SEED/JSEED,IFRST,NEXTN |
461 |
INTEGER JSEED,IFRST,NEXTN |
462 |
INTEGER HVLUE, LVLUE, TESTV |
463 |
C |
464 |
IF (IFRST .EQ. 0) THEN |
465 |
NEXTN = JSEED |
466 |
IFRST = 1 |
467 |
ENDIF |
468 |
C |
469 |
HVLUE = NEXTN / MOBYMP |
470 |
LVLUE = MOD(NEXTN, MOBYMP) |
471 |
TESTV = MPLIER*LVLUE - MOMDMP*HVLUE |
472 |
IF (TESTV .GT. 0) THEN |
473 |
NEXTN = TESTV |
474 |
ELSE |
475 |
NEXTN = TESTV + MODLUS |
476 |
ENDIF |
477 |
RAND = REAL(NEXTN)/REAL(MODLUS) |
478 |
C |
479 |
RETURN |
480 |
END |
481 |
#endif |