Context Navigation

physiq.F @ 135

Last change on this file since 135 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 47.7 KB

Line
1	subroutine physiq(ngrid,nlayer,nq,
2	* firstcall,lastcall,
3	* pday,ptime,ptimestep,
4	* pplev,pplay,pphi,
5	* pu,pv,pt,pq,
6	* pw,
7	* pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
8
9	use radinc_h, only : naerkind
10
11	implicit none
12
13
14	!==================================================================
15	!
16	! Purpose
17	! -------
18	! Central subroutine for all the physics parameterisations in the
19	! universal model. Originally adapted from the Mars LMDZ model.
20	!
21	! The model can be run without or with tracer transport
22	! depending on the value of "tracer" in file "callphys.def".
23	!
24	!
25	! It includes:
26	!
27	! 1. Initialization:
28	! 1.1 Firstcall initializations
29	! 1.2 Initialization for every call to physiq
30	! 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
31	! 2. Compute radiative transfer tendencies
32	! (longwave and shortwave).
33	! 4. Vertical diffusion (turbulent mixing):
34	! 5. Convective adjustment
35	! 6. Condensation and sublimation of gases (currently just CO2).
36	! 7. TRACERS :
37	! 7a. water and water ice
38	! 7c. other schemes for tracer transport (lifting, sedimentation)
39	! 7d. updates (pressure variations, surface budget)
40	! 9. Surface and sub-surface temperature calculations
41	! 10. Write outputs :
42	! - "startfi", "histfi" if it's time
43	! - Saving statistics if "callstats = .true."
44	! - Output any needed variables in "diagfi"
45	! 10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
46	!
47	! arguments
48	! ---------
49	!
50	! input
51	! -----
52	! ecri period (in dynamical timestep) to write output
53	! ngrid Size of the horizontal grid.
54	! All internal loops are performed on that grid.
55	! nlayer Number of vertical layers.
56	! nq Number of advected fields
57	! firstcall True at the first call
58	! lastcall True at the last call
59	! pday Number of days counted from the North. Spring
60	! equinoxe.
61	! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
62	! ptimestep timestep (s)
63	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
64	! pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
65	! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
66	! pu(ngrid,nlayer) u component of the wind (ms-1)
67	! pv(ngrid,nlayer) v component of the wind (ms-1)
68	! pt(ngrid,nlayer) Temperature (K)
69	! pq(ngrid,nlayer,nq) Advected fields
70	! pudyn(ngrid,nlayer) \
71	! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
72	! ptdyn(ngrid,nlayer) / corresponding variables
73	! pqdyn(ngrid,nlayer,nq) /
74	! pw(ngrid,?) vertical velocity
75	!
76	! output
77	! ------
78	!
79	! pdu(ngrid,nlayermx) \
80	! pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
81	! pdt(ngrid,nlayermx) / variables due to physical processes.
82	! pdq(ngrid,nlayermx) /
83	! pdpsrf(ngrid) /
84	! tracerdyn call tracer in dynamical part of GCM ?
85	!
86	!
87	! Authors
88	! -------
89	! Frederic Hourdin 15/10/93
90	! Francois Forget 1994
91	! Christophe Hourdin 02/1997
92	! Subroutine completly rewritten by F.Forget (01/2000)
93	! Water ice clouds: Franck Montmessin (update 06/2003)
94	! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
95	! New correlated-k radiative scheme: R. Wordsworth (2009)
96	! Many specifically Martian subroutines removed: R. Wordsworth (2009)
97	!
98	!==================================================================
99
100
101	c 0. Declarations :
102	c ------------------
103
104	#include "dimensions.h"
105	#include "dimphys.h"
106	#include "comgeomfi.h"
107	#include "surfdat.h"
108	#include "comsoil.h"
109	#include "comdiurn.h"
110	#include "callkeys.h"
111	#include "comcstfi.h"
112	#include "planete.h"
113	#include "comsaison.h"
114	#include "control.h"
115	#include "comg1d.h"
116	#include "tracer.h"
117
118	#include "watercap.h"
119
120	#include "netcdf.inc"
121
122
123
124	c Arguments :
125	c -----------
126
127	c inputs:
128	c -------
129	INTEGER ngrid,nlayer,nq
130	REAL ptimestep
131	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
132	REAL pphi(ngridmx,nlayer)
133	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
134	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
135	REAL pw(ngridmx,nlayer) ! pvervel transmitted by dyn3d
136	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
137	LOGICAL firstcall,lastcall
138
139	REAL pday
140	REAL ptime
141	logical tracerdyn
142
143	c outputs:
144	c --------
145	c physical tendencies
146	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
147	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
148	REAL pdpsrf(ngridmx) ! surface pressure tendency
149
150
151	c Local saved variables:
152	c ----------------------
153	c aerosol (dust or ice) extinction optical depth at reference wavelength
154	c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
155	c aerosol optical properties:
156	c REAL aerosol(ngridmx,nlayermx,naerkind)
157	! this is now internal to callcorrk and hence no longer needed here
158
159	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
160	INTEGER icount ! counter of calls to physiq during the run.
161	REAL tsurf(ngridmx) ! Surface temperature (K)
162	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
163	REAL albedo(ngridmx) ! Surface albedo
164
165	real albedo0(ngridmx) ! Surface albedo
166	save albedo0
167
168	REAL emis(ngridmx) ! Thermal IR surface emissivity
169	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
170	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
171	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
172	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
173	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
174	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
175	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
176
177	SAVE day_ini, icount
178	SAVE tsurf,tsoil
179	SAVE albedo,emis, q2
180	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
181
182	REAL stephan
183	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
184	SAVE stephan
185
186	c Local variables :
187	c -----------------
188
189	! aerosol (dust or ice) extinction optical depth at reference wavelength
190	! for the "naerkind" optically active aerosols:
191	REAL aerosol(ngridmx,nlayermx,naerkind)
192
193	CHARACTER*80 fichier
194	INTEGER l,ig,ierr,igout,iq,i, tapphys
195	! INTEGER iqmin ! Used if iceparty engaged
196
197	REAL fluxsurf_lw(ngridmx) ! incident LW (IR) surface flux (W.m-2)
198	REAL fluxsurf_sw(ngridmx) ! incident SW (solar) surface flux (W.m-2)
199	REAL fluxtop_lw(ngridmx) ! Outgoing LW (IR) flux to space (W.m-2)
200	REAL fluxtop_sw(ngridmx) ! Outgoing SW (solar) flux to space (W.m-2)
201
202	! included by RW for equilibration test
203	real fluxtop_dn(ngridmx)
204	real fluxabs_sw(ngridmx) ! absorbed shortwave radiation
205	real fluxdyn(ngridmx) ! horizontal heat transport by dynamics
206
207	REAL zls ! solar longitude (rad)
208	REAL zday ! date (time since Ls=0, in martian days)
209	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
210	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
211	REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
212
213	c Aerosol effective radius used for radiative transfer (units=meters)
214	REAL reffrad(ngridmx,nlayermx,naerkind)
215
216	c Tendencies due to various processes:
217	REAL dqsurf(ngridmx,nqmx)
218	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
219	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
220	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear areas
221	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear areas
222	REAL zdtsurf(ngridmx) ! (K/s)
223	REAL zdtlscale(ngridmx,nlayermx)
224	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
225	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
226	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
227	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
228	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
229	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
230	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
231	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
232
233	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
234	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
235	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
236	REAL zdqadj(ngridmx,nlayermx,nqmx)
237	REAL zdqc(ngridmx,nlayermx,nqmx)
238	REAL zdqlscale(ngridmx,nlayermx,nqmx)
239	REAL zdqslscale(ngridmx,nqmx)
240	REAL zdqchim(ngridmx,nlayermx,nqmx)
241	REAL zdqschim(ngridmx,nqmx)
242
243	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
244	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
245	REAL zdumolvis(ngridmx,nlayermx)
246	REAL zdvmolvis(ngridmx,nlayermx)
247	real zdqmoldiff(ngridmx,nlayermx,nqmx)
248
249	c Local variable for local intermediate calcul:
250	REAL zflubid(ngridmx)
251	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
252	REAL zdum1(ngridmx,nlayermx)
253	REAL zdum2(ngridmx,nlayermx)
254	REAL ztim1,ztim2,ztim3, z1,z2
255	REAL ztime_fin
256	REAL zdh(ngridmx,nlayermx)
257	INTEGER length
258	PARAMETER (length=100)
259
260	c local variables only used for diagnostic (output in file "diagfi" or "stats")
261	c -----------------------------------------------------------------------------
262	REAL ps(ngridmx), zt(ngridmx,nlayermx)
263	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
264	REAL zq(ngridmx,nlayermx,nqmx)
265	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
266	character*2 str2
267	character*5 str5
268	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
269	real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
270	save ztprevious
271	real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
272	real qtot1,qtot2 ! total aerosol mass
273	integer igmin, lmin
274	logical tdiag
275
276	real co2col(ngridmx) ! CO2 column
277	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
278	REAL zstress(ngrid), cd
279	real hco2(nqmx),tmean, zlocal(nlayermx)
280	real rho(ngridmx,nlayermx) ! density
281	real vmr(ngridmx,nlayermx) ! volume mixing ratio
282
283	REAL time_phys
284
285	! included by RW for kastprof scheme
286	logical kastprof
287
288	! reinstated by RW for diagnostic
289	REAL tau_col(ngrid)
290
291	! included by RW for tidally locked dynamics
292	! logical tlocked
293
294	! included by RW to reduce insanity of code
295	real flux1,flux2,flux3,ts1,ts2,ts3
296
297	! included by RW for temporary comparison
298	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
299
300	! included by RW to compute tracer column densities
301	real qcol(ngridmx,nqmx)
302
303	! included by RW for H2O precipitation
304	REAL zdtrain(ngridmx,nlayermx)
305	REAL zdqrain(ngridmx,nlayermx,nqmx)
306	REAL zdqsrain(ngridmx)
307	REAL zdqssnow(ngridmx)
308	REAL rainout(ngridmx)
309
310	! included by RW for H2O Manabe scheme
311	REAL dtmoist(ngridmx,nlayermx)
312	REAL dqmoist(ngridmx,nlayermx,nqmx)
313
314	real qvap(ngridmx,nlayermx)
315	real dqvaplscale(ngridmx,nlayermx)
316	real dqcldlscale(ngridmx,nlayermx)
317	real cloudfrac(ngridmx,nlayermx)
318	real rneb_man(ngridmx,nlayermx)
319	real rneb_lsc(ngridmx,nlayermx)
320
321	! included by RW to account for surface cooling by evaporation
322	REAL dtsurfh2olat(ngridmx)
323
324	! included by RW to test energy conservation
325	REAL dEtot, dEtots, masse, vabs, dvabs
326
327	! included by RW to test water conservation
328	real h2otot
329
330	! included by RW to allow variations in cp with location
331	REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure
332	REAL rcp3D(ngridmx,nlayermx) ! R / specific heat capacity
333	real cppNI, rcpNI, nnu ! last one just for Seb version
334	REAL zpopskNI(ngridmx,nlayermx)
335
336	! included by RW to make 1D saves not every timestep
337	integer countG1D
338	save countG1D
339
340	c=======================================================================
341
342	kastprof=.false.
343
344	c 1. Initialisation:
345	c -----------------
346
347	c 1.1 Initialisation only at first call
348	c ---------------------------------------
349	IF (firstcall) THEN
350
351	if(ngrid.eq.1)then
352	!saveG1D=day_step*10
353	saveG1D=1
354	countG1D=1
355	endif
356
357	c variables set to 0
358	c ~~~~~~~~~~~~~~~~~~
359	call zerophys(ngrid*nlayer,dtrad)
360	call zerophys(ngrid,fluxrad)
361
362	c initialize tracer names, indexes and properties
363	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
364	tracerdyn=tracer
365	IF (tracer) THEN
366	CALL initracer()
367	ENDIF ! end tracer
368
369
370	c read startfi (initial parameters)
371	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
372
373	CALL phyetat0 ("startfi.nc",0,0,
374	& nsoilmx,nq,
375	& day_ini,time_phys,
376	& tsurf,tsoil,emis,q2,qsurf)
377
378	if (pday.ne.day_ini) then
379	write(,) "PHYSIQ: ERROR: bad synchronization between ",
380	& "physics and dynamics"
381	write(,) "dynamics day: ",pday
382	write(,) "physics day: ",day_ini
383	stop
384	endif
385
386	write (,) 'In physiq day_ini =', day_ini
387
388
389	c Initialize albedo and orbital calculation
390	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
391	CALL surfini(ngrid,qsurf,albedo)
392	! surely for good CO2 ice feedback this must be called every turn?
393	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
394
395
396	if(tlocked)then
397	print*,'Planet is tidally locked at resonance n=',nres
398	print*,'Make sure you have the right rotation rate!!!'
399	endif
400
401
402	c initialize soil
403	c ~~~~~~~~~~~~~~~
404	IF (callsoil) THEN
405	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
406	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
407	ELSE
408	PRINT*,'WARNING! Thermal conduction in the soil turned off'
409	DO ig=1,ngrid
410	capcal(ig)=1.e6
411	fluxgrd(ig)=0.
412	ENDDO
413	ENDIF
414	icount=1
415
416	c Set temperature just above condensation temperature (for Early Mars)
417	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
418	if(nearco2cond) then
419	write(,)' WARNING, starting at Tcond+1K !!!'
420	DO l=1, nlayer
421	DO ig=1,ngrid
422	pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4
423	& -pt(ig,l)) / ptimestep
424	ENDDO
425	ENDDO
426	endif
427
428	if(meanOLR)then
429	! to record global radiative balance
430	call system('rm -f rad_bal.out')
431	endif
432
433	if(water)then
434	! initialise variables for water cycle
435
436	if(ngrid.eq.1)then
437	qsurf(1,igcm_h2o_ice) = 100.0
438	DO l=1, nlayer
439	pq(1,l,igcm_h2o_ice)=0.0!1.0e-2
440	enddo
441	endif
442
443	endif
444	call su_watercycle ! even if we don't have a water cycle, we might
445	! need epsi for the wvp definitions in callcorrk.F
446
447	ENDIF ! (end of "if firstcall")
448
449	c ---------------------------------------------------
450	c 1.2 Initializations done at every physical timestep:
451	c ---------------------------------------------------
452	c
453	IF (ngrid.NE.ngridmx) THEN
454	PRINT*,'STOP in PHYSIQ'
455	PRINT*,'Probleme de dimensions :'
456	PRINT*,'ngrid = ',ngrid
457	PRINT*,'ngridmx = ',ngridmx
458	STOP
459	ENDIF
460
461	c Initialize various variables
462	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
463	call zerophys(ngrid*nlayer, pdv)
464	call zerophys(ngrid*nlayer, pdu)
465	if ( (.not.nearco2cond).and.(.not.firstcall) ) then
466	call zerophys(ngrid*nlayer,pdt)
467	endif ! this was the source of an evil bug...
468	call zerophys(ngridnlayernq, pdq)
469	call zerophys(ngrid, pdpsrf)
470	call zerophys(ngrid, zflubid)
471	call zerophys(ngrid, zdtsurf)
472	call zerophys(ngrid*nq, dqsurf)
473	igout=ngrid/2+1
474
475	zday=pday+ptime ! compute time, in sols (and fraction thereof)
476
477	c Compute Solar Longitude (Ls) :
478	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
479	if (season) then
480	call solarlong(zday,zls)
481	else
482	call solarlong(float(day_ini),zls)
483	end if
484
485	c Compute geopotential between layers
486	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
487
488	DO l=1,nlayer
489	DO ig=1,ngrid
490	zzlay(ig,l)=pphi(ig,l)/g
491	ENDDO
492	ENDDO
493	DO ig=1,ngrid
494	zzlev(ig,1)=0.
495	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
496	ENDDO
497	DO l=2,nlayer
498	DO ig=1,ngrid
499	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
500	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
501	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
502	ENDDO
503	ENDDO
504
505	! Potential temperature calculation not the same in physiq and dynamic...
506
507	c Compute potential temperature
508	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
509	if(nonideal)then
510
511	DO l=1,nlayer
512	DO ig=1,ngrid
513	call calc_cpp3d(cppNI,rcpNI,pt(ig,l),pplay(ig,l))
514	cpp3D(ig,l) = cppNI
515	rcp3D(ig,l) = rcpNI
516	ENDDO
517	ENDDO
518
519	DO l=1,nlayer
520	DO ig=1,ngrid
521
522	! nnu=0.35
523	! zh(ig,l) = (pt(ig,l)**nnu +
524	! & nnu(460nnu)(8.314/1000.)*
525	! & log(pplev(ig,1)/pplay(ig,l)) )**(-nnu) ! Sebastian's version
526	! zpopskNI(ig,l) = (pplay(ig,l)/pplev(ig,1))**rcp3D(ig,l)
527
528	zh(ig,l) = pt(ig,l)*1.0024 - 0.659
529	zpopskNI(ig,l) = pt(ig,l)/zh(ig,l)
530	! we're only after zpopskNI here, not zh
531	! zh is calculated seperately before both vdifc and convadj
532	ENDDO
533	ENDDO
534
535	endif
536	! else
537
538	DO l=1,nlayer
539	DO ig=1,ngrid
540	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
541	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
542	ENDDO
543	ENDDO
544
545	! endif
546
547
548	c-----------------------------------------------------------------------
549	c 2. Compute radiative tendencies :
550	c-----------------------------------------------------------------------
551
552	IF (callrad) THEN
553	IF( MOD(icount-1,iradia).EQ.0) THEN
554
555	c Local stellar zenith angle
556	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
557	CALL orbite(zls,dist_sol,declin)
558
559	IF (tlocked) THEN
560	ztim1=SIN(declin)
561	ztim2=COS(declin)COS(2.pi(zday/year_day) - zlsnres)
562	ztim3=-COS(declin)SIN(2.pi(zday/year_day) - zlsnres)
563
564	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
565	s ztim1,ztim2,ztim3,mu0,fract)
566
567	ELSEIF (diurnal) THEN
568	ztim1=SIN(declin)
569	ztim2=COS(declin)COS(2.pi*(zday-.5))
570	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
571
572	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
573	s ztim1,ztim2,ztim3,mu0,fract)
574
575	ELSE
576
577	CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
578	! WARNING: this function appears not to work in 1D
579
580	ENDIF
581
582	c Call main radiative transfer scheme
583	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
584
585	c Radiative transfer
586	c ------------------
587
588	if (corrk) then
589	if(kastprof)then
590	call kastprof_fn(tsurf,pt,pplay,pplev)
591	print*,'pt',pt
592	print*,'pplay',pplay
593	endif
594
595	CALL callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,
596	& albedo,emis,mu0,pplev,pplay,pt,
597	& tsurf,fract,dist_sol,igout,aerosol,cpp3D,
598	& zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
599	& fluxtop_sw,fluxtop_dn,reffrad,tau_col,ptime,pday,
600	& firstcall,lastcall)
601
602
603
604	c Radiative flux from the sky absorbed by the surface (W.m-2)
605	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
606	DO ig=1,ngrid
607	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
608	& +fluxsurf_sw(ig)*(1.-albedo(ig))
609	ENDDO
610
611	c Net atmospheric radiative heating rate (K.s-1)
612	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
613
614	DO l=1,nlayer
615	DO ig=1,ngrid
616	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
617	ENDDO
618	ENDDO
619
620	else
621
622	c Atmosphere has no radiative effect
623	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
624	DO ig=1,ngrid
625	fluxtop_dn(ig) = fract(ig)mu0(ig)Fat1AU/dist_sol**2
626	fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
627	fluxtop_sw(ig) = fluxtop_dn(ig)*albedo(ig)
628	fluxtop_lw(ig) = emis(ig)stephantsurf(ig)**4
629	ENDDO ! radiation skips the atmosphere entirely
630
631
632	DO l=1,nlayer
633	DO ig=1,ngrid
634	dtrad(ig,l)=0.0
635	ENDDO
636	ENDDO ! hence no atmospheric radiative heating
637
638	endif ! if corrk
639
640	ENDIF ! of if(mod(icount-1,iradia).eq.0)
641
642	! Transformation of the radiative tendencies:
643	! -------------------------------------------
644
645	! Net radiative surface flux (W.m-2)
646	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
647
648	DO ig=1,ngrid
649	zplanck(ig)=tsurf(ig)*tsurf(ig)
650	zplanck(ig)=emis(ig)*
651	& stephanzplanck(ig)zplanck(ig)
652	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
653	ENDDO
654
655	DO l=1,nlayer
656	DO ig=1,ngrid
657	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
658	ENDDO
659	ENDDO
660
661	!-------------------------
662	! test energy conservation
663	if(enertest)then
664	dEtot=0.0
665	dEtots=0.0
666	DO ig = 1, ngrid
667	DO l = 1, nlayer
668	masse = (pplev(ig,l) - pplev(ig,l+1))/g
669	dEtot = dEtot + cppmassedtrad(ig,l)*area(ig)
670	!print,'l=',l,'dEtot=',cppmassedtrad(ig,l)area(ig)
671	ENDDO
672	dEtots = dEtots + fluxrad(ig)*area(ig)
673	ENDDO
674	dEtot=dEtot/totarea
675	dEtots=dEtots/totarea
676	print*,'In corrk atmospheric energy change=',dEtot,' W m-2'
677	print*,'In corrk surface energy change=',dEtots,' W m-2'
678	endif
679	!-------------------------
680
681	ENDIF ! of IF (callrad)
682
683
684	!-----------------------------------------------------------------------
685	! 4. Vertical diffusion (turbulent mixing):
686	! -----------------------------------------
687
688	IF (calldifv) THEN
689
690	DO ig=1,ngrid
691	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
692	ENDDO
693
694	CALL zerophys(ngrid*nlayer,zdum1)
695	CALL zerophys(ngrid*nlayer,zdum2)
696	do l=1,nlayer
697	do ig=1,ngrid
698	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
699	enddo
700	enddo
701
702	c Calling vdif (Martian version WITH CO2 condensation)
703	CALL vdifc(ngrid,nlayer,nq,zpopsk,
704	& ptimestep,capcal,lwrite,
705	& pplay,pplev,zzlay,zzlev,z0,
706	& pu,pv,zh,pq,tsurf,emis,qsurf,
707	& zdum1,zdum2,zdh,pdq,zflubid,
708	& zdudif,zdvdif,zdhdif,zdtsdif,q2,
709	& zdqdif,zdqsdif)
710
711	DO l=1,nlayer
712	DO ig=1,ngrid
713	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
714	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
715	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
716
717	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
718	ENDDO
719	ENDDO
720
721	DO ig=1,ngrid
722	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
723	ENDDO
724
725	if (tracer) then
726	DO iq=1, nq
727	DO l=1,nlayer
728	DO ig=1,ngrid
729	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
730	ENDDO
731	ENDDO
732	ENDDO
733	DO iq=1, nq
734	DO ig=1,ngrid
735	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
736	ENDDO
737	ENDDO
738
739	end if ! of if (tracer)
740
741	!print*,'surface temp loss:'
742	!print*, dtsurfh2olat
743
744	!-------------------------
745	! test energy conservation
746	if(enertest)then
747	dEtot=0.0
748	dEtots=0.0
749	DO ig = 1, ngrid
750	DO l = 1, nlayer
751	masse = (pplev(ig,l) - pplev(ig,l+1))/g
752
753	dEtot = dEtot + cppmassezdtdif(ig,l)*area(ig)
754
755	vabs = sqrt(pdu(ig,l)2 + pdv(ig,l)2)
756	dvabs = sqrt(zdudif(ig,l)2 + zdvdif(ig,l)2)
757	dEtot = dEtot + massevabsdvabs*area(ig)
758	ENDDO
759	dEtots = dEtots + capcal(ig)zdtsdif(ig)area(ig)
760	ENDDO
761	dEtot=dEtot/totarea
762	dEtots=dEtots/totarea
763	print*,'In difv atmospheric energy change=',dEtot,' W m-2'
764	print*,'In difv surface energy change=',dEtots,' W m-2'
765	endif
766	!-------------------------
767
768	ELSE
769
770	DO ig=1,ngrid
771	zdtsurf(ig)=zdtsurf(ig)+
772	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
773	ENDDO
774
775	ENDIF ! of IF (calldifv)
776
777
778	!-----------------------------------------------------------------------
779	! 5. Dry convective adjustment:
780	! -----------------------------
781
782	IF(calladj) THEN
783
784	DO l=1,nlayer
785	DO ig=1,ngrid
786	if(nonideal)then
787	zdh(ig,l)=pdt(ig,l)/zpopskNI(ig,l)
788	else
789	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
790	endif
791	ENDDO
792	ENDDO
793	CALL zerophys(ngrid*nlayer,zduadj)
794	CALL zerophys(ngrid*nlayer,zdvadj)
795	CALL zerophys(ngrid*nlayer,zdhadj)
796
797	if(nonideal)then
798	CALL convadj(ngrid,nlayer,nq,ptimestep,
799	& pplay,pplev,zpopskNI,
800	& pu,pv,zh,pq,
801	& pdu,pdv,zdh,pdq,
802	& zduadj,zdvadj,zdhadj,
803	& zdqadj)
804	else
805	CALL convadj(ngrid,nlayer,nq,ptimestep,
806	& pplay,pplev,zpopsk,
807	& pu,pv,zh,pq,
808	& pdu,pdv,zdh,pdq,
809	& zduadj,zdvadj,zdhadj,
810	& zdqadj)
811	endif
812
813
814	DO l=1,nlayer
815	DO ig=1,ngrid
816	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
817	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
818	if(nonideal)then
819	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopskNI(ig,l)
820	zdtadj(ig,l)=zdhadj(ig,l)*zpopskNI(ig,l) ! for diagnostic only
821	else
822	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
823	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
824	endif
825	ENDDO
826	ENDDO
827
828	if(tracer) then
829	DO iq=1, nq
830	DO l=1,nlayer
831	DO ig=1,ngrid
832	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
833	ENDDO
834	ENDDO
835	ENDDO
836	end if
837
838	!-------------------------
839	! test energy conservation
840	if(enertest)then
841	dEtot=0.0
842	DO ig = 1, ngrid
843	DO l = 1, nlayer
844	masse = (pplev(ig,l) - pplev(ig,l+1))/g
845	dEtot = dEtot + cppmassezdtadj(ig,l)*area(ig)
846	ENDDO
847	ENDDO
848	dEtot=dEtot/totarea
849	print*,'In convadj atmospheric energy change=',dEtot,' W m-2'
850	endif
851	!-------------------------
852
853	ENDIF ! of IF(calladj)
854
855	!-----------------------------------------------------------------------
856	! 6. Carbon dioxide condensation-sublimation:
857	! -------------------------------------------
858
859	IF (co2cond) THEN
860	if (.not.tracer) then
861	print*,'We need a CO2 ice tracer to condense CO2'
862	call abort
863	endif
864
865	call condens_co2cloud(ngrid,nlayer,nq,ptimestep,
866	& capcal,pplay,pplev,tsurf,pt,
867	& pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
868	& qsurf(1,igcm_co2_ice),albedo,emis,
869	& zdtc,zdtsurfc,pdpsrf,zduc,zdvc,
870	& zdqc,reffrad,cpp3D)
871
872	DO l=1,nlayer
873	DO ig=1,ngrid
874	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
875	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
876	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
877	ENDDO
878	ENDDO
879	DO ig=1,ngrid
880	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
881	ENDDO
882
883	! IF (tracer) THEN
884	DO iq=1, nq ! should use new notation here !
885	DO l=1,nlayer
886	DO ig=1,ngrid
887	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
888	ENDDO
889	ENDDO
890	ENDDO
891	!NB: we do not add tendency on surface co2ice ;
892	! since qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
893	! ENDIF ! of IF (tracer)
894
895
896	!-------------------------
897	! test energy conservation
898	if(enertest)then
899	dEtot=0.0
900	DO ig = 1, ngrid
901	DO l = 1, nlayer
902	masse = (pplev(ig,l) - pplev(ig,l+1))/g
903	dEtot = dEtot + cppmassezdtc(ig,l)
904	ENDDO
905	dEtot = dEtot + capcal(ig)*zdtsurfc(ig)
906	ENDDO
907	print*,'In co2cloud energy change=',dEtot,' W m-2'
908	endif ! we're missing the latent heat...
909	!-------------------------
910
911	ENDIF ! of IF (co2cond)
912
913
914
915	c-----------------------------------------------------------------------
916	c 7. Specific parameterizations for tracers
917	c -----------------------------------------
918
919	if (tracer) then
920
921	c 7a. Water and ice
922	c ---------------
923
924	c ---------------------------------------
925	c Water ice condensation in the atmosphere
926	c ----------------------------------------
927	IF (water) THEN
928
929	if(watercond)then
930
931	if(1.eq.2)then
932	call moistadj(pt,pq,pplev,pplay,dtmoist,dqmoist,
933	& ptimestep,rneb_man)
934
935	DO l=1,nlayer
936	DO ig=1,ngrid
937	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
938	& dqmoist(ig,l,igcm_h2o_vap)
939	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
940	& dqmoist(ig,l,igcm_h2o_ice)
941	pdt(ig,l)=pdt(ig,l)+dtmoist(ig,l)
942	ENDDO
943	ENDDO
944	endif ! moistadj turned off for now
945
946	DO l=1,nlayer
947	DO ig=1,ngrid
948	qvap(ig,l)=pq(ig,l,igcm_h2o_vap)
949	ENDDO
950	ENDDO
951	call largescale2(ptimestep,pplev,pplay,pt,qvap,
952	s zdtlscale, dqvaplscale,dqcldlscale,rneb_lsc)
953	DO l=1,nlayer
954	DO ig=1,ngrid
955	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
956	& dqvaplscale(ig,l)
957	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
958	& dqcldlscale(ig,l)
959	! pdt(ig,l)=pdt(ig,l)+zdtlscale(ig,l)
960	! for the moment this is removed, as it dont conserve energy
961	ENDDO
962	ENDDO
963
964	! compute cloud fraction
965	DO l = 1, nlayer
966	DO ig = 1,ngrid
967	cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
968	ENDDO
969	ENDDO
970
971	endif ! of IF (watercondense)
972
973	if(waterrain)then
974
975	call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,
976	& zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
977
978	DO l=1,nlayer
979	DO ig=1,ngrid
980	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
981	& zdqrain(ig,l,igcm_h2o_vap)
982	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
983	& zdqrain(ig,l,igcm_h2o_ice)
984	!pdt(ig,l)=pdt(ig,l)+zdtrain(ig,l)
985	ENDDO
986	ENDDO
987
988	DO ig=1,ngrid
989	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
990	& zdqsrain(ig)+zdqssnow(ig)
991	rainout(ig)=zdqsrain(ig)+zdqssnow(ig) ! diagnostic
992	!rainout(ig)=zdqssnow(ig)
993	ENDDO
994
995	end if ! of IF (waterrain)
996
997	end if ! of IF (water)
998
999
1000	c 7c. Aerosol particles
1001	c -------------------
1002
1003	c -------------
1004	c Sedimentation
1005	c -------------
1006	IF (sedimentation) THEN
1007	call zerophys(ngridnlayernq, zdqsed)
1008	call zerophys(ngrid*nq, zdqssed)
1009
1010	call callsedim(ngrid,nlayer, ptimestep,
1011	& pplev,zzlev, pt,
1012	& pq, pdq, zdqsed, zdqssed,nq)
1013
1014	DO iq=2, nq ! should be updated
1015	DO l=1,nlayer
1016	DO ig=1,ngrid
1017	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1018	ENDDO
1019	ENDDO
1020	ENDDO
1021	DO iq=2, nq
1022	DO ig=1,ngrid
1023	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1024	ENDDO
1025	ENDDO
1026	END IF ! of IF (sedimentation)
1027
1028
1029	c 7d. Updates
1030	c ---------
1031
1032	DO iq=1, nq
1033	DO ig=1,ngrid
1034
1035	c ---------------------------------
1036	c Updating tracer budget on surface
1037	c ---------------------------------
1038	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1039
1040	ENDDO ! (ig)
1041	ENDDO ! (iq)
1042
1043	endif ! of if (tracer)
1044
1045
1046	!-----------------------------------------------------------------------
1047	! 9. Surface and sub-surface soil temperature
1048	!-----------------------------------------------------------------------
1049	!
1050	!
1051	! 9.1 Increment Surface temperature:
1052	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1053
1054	DO ig=1,ngrid
1055	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1056	ENDDO
1057
1058	!
1059	! 9.2 Compute soil temperatures and subsurface heat flux:
1060	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1061	IF (callsoil) THEN
1062	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1063	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
1064	ENDIF
1065
1066	!-------------------------
1067	! test energy conservation
1068	if(enertest)then
1069	dEtots=0.0
1070	DO ig = 1, ngrid
1071	dEtots = dEtots + capcal(ig)zdtsurf(ig)area(ig)
1072	ENDDO
1073	dEtots=dEtots/totarea
1074	print*,'Surface energy change=',dEtots,' W m-2'
1075	endif
1076	!-------------------------
1077
1078
1079
1080
1081	!-----------------------------------------------------------------------
1082	! 10. Perform diagnostics and write output files
1083	!-----------------------------------------------------------------------
1084
1085	! ---------------------------------------------------------
1086	! Check the energy balance of the simulation during the run
1087	! ---------------------------------------------------------
1088	flux1 = 0
1089	flux2 = 0
1090	flux3 = 0
1091	do ig=1,ngrid
1092
1093	flux1 = flux1 +
1094	& area(ig)*(fluxtop_dn(ig) - fluxtop_sw(ig))
1095	flux2 = flux2 + area(ig)*fluxtop_lw(ig)
1096	flux3 = flux3 + area(ig)*fluxtop_dn(ig)
1097	fluxabs_sw(ig)=fluxtop_dn(ig) - fluxtop_sw(ig)
1098
1099	end do
1100	print*,'Incident solar flux, absorbed solar flux, OLR (W/m2)'
1101	print*, flux3/totarea, ' ', flux1/totarea ,
1102	& ' = ', flux2/totarea
1103
1104	if(meanOLR)then
1105	! to record global radiative balance
1106	open(92,file="rad_bal.out",form='formatted',access='append')
1107	write(92,*) zday,flux3/totarea,flux1/totarea,flux2/totarea
1108	close(92)
1109	endif
1110
1111	if(kastprof)then ! for f(p,T) computation (not needed in universal model)
1112	open(15,file='ene_bal.out',form='formatted',access='append')
1113	write(15,*) flux1
1114	close(15)
1115	stop
1116	endif
1117
1118	! Surface temperature information
1119	ts1 = 0
1120	ts2 = 999
1121	ts3 = 0
1122	do ig=1,ngrid
1123	ts1 = ts1 + area(ig)*tsurf(ig)
1124	ts2 = min(ts2,tsurf(ig))
1125	ts3 = max(ts3,tsurf(ig))
1126	end do
1127	print*,'Mean Tsurf =',ts1/totarea , ' Min Tsurf=',ts2,
1128	& ' Max Tsurf =',ts3
1129
1130	! -------------------------------
1131	! Dynamical fields incrementation
1132	! -------------------------------
1133	! (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1134	! temperature, zonal and meridional wind
1135	DO l=1,nlayer
1136	DO ig=1,ngrid
1137	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1138	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1139	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1140
1141	! diagnostic
1142	zdtdyn(ig,l)=ztprevious(ig,l)-pt(ig,l)
1143	ztprevious(ig,l)=zt(ig,l)
1144	ENDDO
1145	ENDDO
1146	if(firstcall)call zerophys(zdtdyn)
1147
1148	! dynamical heating diagnostic
1149	DO ig=1,ngrid
1150	fluxdyn(ig)=0.
1151	DO l=1,nlayer
1152	if(nonideal)then
1153	fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)
1154	& (pplev(ig,l)-pplev(ig,l+1))cpp3D(ig,l)/g
1155	else
1156	fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)
1157	& (pplev(ig,l)-pplev(ig,l+1))cpp/g
1158	endif
1159	ENDDO
1160	ENDDO
1161
1162
1163	! tracers
1164	DO iq=1, nq
1165	DO l=1,nlayer
1166	DO ig=1,ngrid
1167	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1168	ENDDO
1169	ENDDO
1170	ENDDO
1171
1172	! surface pressure
1173	DO ig=1,ngrid
1174	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1175	ENDDO
1176
1177	! pressure
1178	DO l=1,nlayer
1179	DO ig=1,ngrid
1180	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1181	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1182	ENDDO
1183	ENDDO
1184
1185
1186	! ---------------------------------------------------------
1187	! Compute column amounts (kg m-2) if tracers are enabled
1188	! ---------------------------------------------------------
1189	if(tracer)then
1190	call zerophys(ngrid*nq,qcol)
1191	do iq=1,nq
1192	DO ig=1,ngrid
1193	DO l=1,nlayer
1194	qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) *
1195	& (pplev(ig,l) - pplev(ig,l+1)) / g
1196	enddo
1197	enddo
1198	enddo
1199	endif
1200
1201
1202	! ---------------------------------------------------------
1203	! Test for water conservation if water is enabled
1204	! ---------------------------------------------------------
1205	if(water)then
1206
1207	h2otot = 0
1208	do ig=1,ngrid
1209	h2otot = h2otot + area(ig)*
1210	& (qcol(ig,igcm_h2o_ice)+qcol(ig,igcm_h2o_vap)
1211	& +qsurf(ig,igcm_h2o_ice)+qsurf(ig,igcm_h2o_vap))
1212	end do
1213	print*,'Total water amount (kg): ',h2otot!/totarea
1214	endif
1215
1216
1217	! TEMPORARY: save water cycle diagnostics
1218	! IF((ngrid.eq.1).and.lastcall) THEN
1219	! open(13,file='Ts.out')
1220	! open(14,file='T.out')
1221	! open(15,file='p.out')
1222	! open(16,file='vap.out')
1223	! open(17,file='alt.out')
1224	! ig=1
1225	! write(13,*) tsurf(ig)
1226	! DO l=1,nlayer
1227	! write(14,*) pt(ig,l)
1228	! write(15,*) pplay(ig,l)
1229	! write(16,*) pq(ig,l,igcm_h2o_vap)
1230	! write(17,*) zzlay(ig,l)
1231	! enddo
1232	! close(13)
1233	! close(14)
1234	! close(15)
1235	! close(16)
1236	! close(17)
1237	! ENDIF
1238
1239
1240	! TEMPORARY: save data for LBL
1241	! IF((ngrid.eq.1).and.lastcall) THEN
1242	! open(14,file='T.dat')
1243	! open(15,file='p.dat')
1244	! ig=1
1245	! write(14,*) nlayer+1
1246	! write(15,*) nlayer+1
1247	! write(14,*) tsurf(ig)
1248	! write(15,*) pplev(ig,1)
1249	! DO l=1,nlayer
1250	! write(14,*) pt(ig,l)
1251	! write(15,*) pplay(ig,l)
1252	! enddo
1253	! close(14)
1254	! close(15)
1255	! ENDIF
1256
1257
1258	IF (ngrid.NE.1) THEN
1259	print,'Ls =',zls180./pi
1260
1261	! -------------------------------------------------------------------
1262	! Writing NetCDF file "RESTARTFI" at the end of the run
1263	! -------------------------------------------------------------------
1264	! Note: 'restartfi' is stored just before dynamics are stored
1265	! in 'restart'. Between now and the writting of 'restart',
1266	! there will have been the itau=itau+1 instruction and
1267	! a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1268	! thus we store for time=time+dtvr
1269
1270	IF(lastcall) THEN
1271	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1272
1273	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1274	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1275	. ptimestep,pday,
1276	. ztime_fin,tsurf,tsoil,emis,q2,qsurf,
1277	. area,albedodat,inertiedat,zmea,zstd,zsig,
1278	. zgam,zthe)
1279	ENDIF
1280
1281	! -----------------------------------------------------------------
1282	! Saving statistics :
1283	! -----------------------------------------------------------------
1284	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
1285	! which can later be used to make the statistic files of the run:
1286	! "stats") only possible in 3D runs !
1287
1288
1289	IF (callstats) THEN
1290
1291	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1292	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1293	call wstats(ngrid,"fluxsurf_lw",
1294	. "Thermal IR radiative flux to surface","W.m-2",2,
1295	. fluxsurf_lw)
1296	call wstats(ngrid,"fluxsurf_sw",
1297	. "Solar radiative flux to surface","W.m-2",2,
1298	. fluxsurf_sw_tot)
1299	call wstats(ngrid,"fluxtop_lw",
1300	. "Thermal IR radiative flux to space","W.m-2",2,
1301	. fluxtop_lw)
1302	call wstats(ngrid,"fluxtop_sw",
1303	. "Solar radiative flux to space","W.m-2",2,
1304	. fluxtop_sw_tot)
1305
1306	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1307	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1308	call wstats(ngrid,"v","Meridional (North-South) wind",
1309	. "m.s-1",3,zv)
1310	call wstats(ngrid,"w","Vertical (down-up) wind",
1311	. "m.s-1",3,pw)
1312	call wstats(ngrid,"q2",
1313	. "Boundary layer eddy kinetic energy","m2.s-2",3,q2)
1314
1315	if (tracer) then
1316	if (water) then
1317	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1318	& *mugaz/mmol(igcm_h2o_vap)
1319	call wstats(ngrid,"vmr_h2ovapor",
1320	. "H2O vapor volume mixing ratio","mol/mol",
1321	. 3,vmr)
1322	endif ! of if (water)
1323
1324	endif !tracer
1325
1326	IF(lastcall) THEN
1327	write (,) "Writing stats..."
1328	call mkstats(ierr)
1329	ENDIF
1330	ENDIF !if callstats
1331
1332
1333
1334	! ----------------------------------------------------------------------
1335	! output in netcdf file "DIAGFI", containing any variable for diagnostic
1336	! (output with period "ecritphy", set in "run.def")
1337	! ----------------------------------------------------------------------
1338	! WRITEDIAGFI can ALSO be called from any other subroutine
1339	! for any variable!
1340
1341	! is it not preferable to keep all the calls in one place?
1342
1343	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1344	& tsurf)
1345	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1346	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1347	! call WRITEDIAGFI(ngrid,"tau","tau"," ",2,tau)
1348	! call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1349	! call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1350	! call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1351	! call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1352	! call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1353	! call WRITEDIAGFI(ngridm,'Teta','T potentielle','K',3,zh)
1354	! call WRITEDIAGFI(ngridm,'Pression','Pression','Pa',3,pplay)
1355
1356	! Total energy balance diagnostics
1357	call WRITEDIAGFI(ngrid,"ISR","incoming stellar rad.","W m-2",
1358	& 2,fluxtop_dn)
1359	call WRITEDIAGFI(ngrid,"ASR","absorbed stellar rad.","W m-2",
1360	& 2,fluxabs_sw)
1361	call WRITEDIAGFI(ngrid,"OLR","outgoing longwave rad.","W m-2",
1362	& 2,fluxtop_lw)
1363	call WRITEDIAGFI(ngrid,"DYN","dynamical heat input","W m-2",
1364	& 2,fluxdyn)
1365
1366	! Temporary inclusions for heating diagnostics
1367	! call WRITEDIAGFI(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
1368	! call WRITEDIAGFI(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
1369	! call WRITEDIAGFI(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
1370
1371	! Output tracers
1372	if (tracer) then
1373	do iq=1,nq
1374	! call WRITEDIAGFI(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
1375	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_surf',
1376	& trim(noms(iq))//'_surf',
1377	& 'kg/kg',2,qsurf(1,iq) )
1378
1379	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_col',
1380	& trim(noms(iq))//'_col',
1381	& 'kg m^-2',2,qcol(1,iq) )
1382
1383
1384
1385	if(waterrain)then
1386	CALL WRITEDIAGFI(ngridmx,"rain",
1387	& "total precipitation","kg m-2 s-1",2,rainout)
1388	endif
1389
1390	! call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_aero',
1391	! & trim(noms(iq))//'_aero',
1392	! & 'kg/kg',3,aerosol(1,1,iq))
1393
1394	call WRITEDIAGFI(ngridmx,"tau_col",
1395	& "Total aerosol optical depth","[]",2,tau_col)
1396
1397	! call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_reff',
1398	! & trim(noms(iq))//'_reff',
1399	! & 'm',3,reffrad(1,1,iq))
1400
1401	enddo
1402	endif
1403
1404	ELSE ! if(ngrid.eq.1)
1405
1406
1407	! ----------------------------------------------------------------------
1408	! Output in grads file "g1d" (ONLY when using testphys1d)
1409	! ----------------------------------------------------------------------
1410
1411	if(countG1D.eq.saveG1D)then
1412
1413	! call WRITEDIAGFI(ngrid,"tsurf","Surface temperature",
1414	! & "K",0,tsurf)
1415	CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
1416	CALL writeg1d(ngrid,1,ps,'ps','Pa')
1417	CALL writeg1d(ngrid,nlayer,zt,'T','K')
1418	CALL writeg1d(ngrid,nlayer,pq(1,1,1),'q','kg/kg')
1419	CALL writeg1d(ngrid,nlayer,aerosol,'aerosol','SI')
1420	CALL writeg1d(ngrid,nlayer,reffrad,'reffrad','SI')
1421	CALL writeg1d(ngrid,nlayer,zdtlw,'dtlw','SI')
1422	CALL writeg1d(ngrid,nlayer,zdtsw,'dtsw','SI')
1423	CALL writeg1d(ngrid,nlayer,zdtdyn,'dtdyn','SI')
1424
1425	! radiation balance (optional)
1426	CALL writeg1d(ngrid,1,flux3/totarea,'ISR','W m-2')
1427	CALL writeg1d(ngrid,1,flux1/totarea,'ASR','W m-2')
1428	CALL writeg1d(ngrid,1,flux2/totarea,'OLR','W m-2')
1429
1430	if(tracer) then
1431	do iq=1,nq
1432	CALL writeg1d(ngrid,1,qsurf(1,iq),
1433	& trim(noms(iq))//'_s','kg m^-2')
1434	CALL writeg1d(ngrid,1,qcol(1,iq),
1435	& trim(noms(iq))//'_c','kg m^-2')
1436	CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
1437	end do
1438
1439	if(waterrain)then
1440	CALL writeg1d(ngrid,1,rainout,
1441	& 'rainfall','kg m-2 s-1')
1442	endif
1443	end if
1444
1445	countG1D=1
1446	else
1447	countG1D=countG1D+1
1448	endif ! if time to save
1449
1450	END IF ! if(ngrid.ne.1)
1451
1452	icount=icount+1
1453
1454	RETURN
1455	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format