Context Navigation

physiq.F @ 260

Last change on this file since 260 was 256, checked in by acolaitis, 14 years ago

--- AC 03/08/2011 ---

M 255 libf/phymars/physiq.F
<> Modified to interface new SL parametrization

M 255 libf/phymars/vdif_cd.F
<> New SL parametrization based on a bulk Richardson Monin-Obukhov theory formulation

Stability functions are taken from D.E. England et al. (95)
Similarity functions coefficients based on Dyer and Hicks (70).
Includes thermal roughness length computation, heat and momentum drag coefficient computation
Can be used to output hydrodynamic-related SL quantities (bulk Richardson, turbulent Prandtl number estimation, Reynolds number)

M 255 libf/phymars/thermcell_dqupdown.F90
<> Minor modification to suit picky compilers

M 255 libf/phymars/vdifc.F
<> Now takes into account sub-grid gustiness, evaluated from thermals activity (it's proxy being the maximum vertical velocity)

M 255 libf/phymars/meso_inc/meso_inc_les.F
<> Minor modification : u* is now taken from the new vdifc and not recomputed from a simple law

File size: 60.2 KB

Line
1	SUBROUTINE physiq(
2	$ ngrid,nlayer,nq
3	$ ,firstcall,lastcall
4	$ ,pday,ptime,ptimestep
5	$ ,pplev,pplay,pphi
6	$ ,pu,pv,pt,pq
7	$ ,pw
8	$ ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn
9	#ifdef MESOSCALE
10	#include "meso_inc/meso_inc_invar.F"
11	#endif
12	$ )
13
14	IMPLICIT NONE
15	c=======================================================================
16	c
17	c subject:
18	c --------
19	c
20	c Organisation of the physical parametrisations of the LMD
21	c martian atmospheric general circulation model.
22	c
23	c The GCM can be run without or with tracer transport
24	c depending on the value of Logical "tracer" in file "callphys.def"
25	c Tracers may be water vapor, ice OR chemical species OR dust particles
26	c
27	c SEE comments in initracer.F about numbering of tracer species...
28	c
29	c It includes:
30	c
31	c 1. Initialization:
32	c 1.1 First call initializations
33	c 1.2 Initialization for every call to physiq
34	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
35	c 2. Compute radiative transfer tendencies
36	c (longwave and shortwave) for CO2 and aerosols.
37	c 3. Gravity wave and subgrid scale topography drag :
38	c 4. Vertical diffusion (turbulent mixing):
39	c 5. Convective adjustment
40	c 6. Condensation and sublimation of carbon dioxide.
41	c 7. TRACERS :
42	c 7a. water and water ice
43	c 7b. call for photochemistry when tracers are chemical species
44	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
45	c 7d. updates (CO2 pressure variations, surface budget)
46	c 8. Contribution to tendencies due to thermosphere
47	c 9. Surface and sub-surface temperature calculations
48	c 10. Write outputs :
49	c - "startfi", "histfi" (if it's time)
50	c - Saving statistics (if "callstats = .true.")
51	c - Dumping eof (if "calleofdump = .true.")
52	c - Output any needed variables in "diagfi"
53	c 11. Diagnostic: mass conservation of tracers
54	c
55	c author:
56	c -------
57	c Frederic Hourdin 15/10/93
58	c Francois Forget 1994
59	c Christophe Hourdin 02/1997
60	c Subroutine completly rewritten by F.Forget (01/2000)
61	c Introduction of the photochemical module: S. Lebonnois (11/2002)
62	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
63	c Water ice clouds: Franck Montmessin (update 06/2003)
64	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
65	c Nb: See callradite.F for more information.
66	c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
67	c
68	c arguments:
69	c ----------
70	c
71	c input:
72	c ------
73	c ecri period (in dynamical timestep) to write output
74	c ngrid Size of the horizontal grid.
75	c All internal loops are performed on that grid.
76	c nlayer Number of vertical layers.
77	c nq Number of advected fields
78	c firstcall True at the first call
79	c lastcall True at the last call
80	c pday Number of days counted from the North. Spring
81	c equinoxe.
82	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
83	c ptimestep timestep (s)
84	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
85	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
86	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
87	c pu(ngrid,nlayer) u component of the wind (ms-1)
88	c pv(ngrid,nlayer) v component of the wind (ms-1)
89	c pt(ngrid,nlayer) Temperature (K)
90	c pq(ngrid,nlayer,nq) Advected fields
91	c pudyn(ngrid,nlayer) \
92	c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
93	c ptdyn(ngrid,nlayer) / corresponding variables
94	c pqdyn(ngrid,nlayer,nq) /
95	c pw(ngrid,?) vertical velocity
96	c
97	c output:
98	c -------
99	c
100	c pdu(ngrid,nlayermx) \
101	c pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
102	c pdt(ngrid,nlayermx) / variables due to physical processes.
103	c pdq(ngrid,nlayermx,nqmx) /
104	c pdpsrf(ngrid) /
105	c tracerdyn call tracer in dynamical part of GCM ?
106
107	c
108	c=======================================================================
109	c
110	c 0. Declarations :
111	c ------------------
112
113	#include "dimensions.h"
114	#include "dimphys.h"
115	#include "comgeomfi.h"
116	#include "surfdat.h"
117	#include "comsoil.h"
118	#include "comdiurn.h"
119	#include "callkeys.h"
120	#include "comcstfi.h"
121	#include "planete.h"
122	#include "comsaison.h"
123	#include "control.h"
124	#include "dimradmars.h"
125	#include "comg1d.h"
126	#include "tracer.h"
127	#include "nlteparams.h"
128
129	#include "chimiedata.h"
130	#include "watercap.h"
131	#include "param.h"
132	#include "param_v3.h"
133	#include "conc.h"
134
135	#include "netcdf.inc"
136
137	#include "slope.h"
138
139	#ifdef MESOSCALE
140	#include "wrf_output_2d.h"
141	#include "wrf_output_3d.h"
142	#include "advtrac.h" !!! this is necessary for tracers (in dyn3d)
143	#include "meso_inc/meso_inc_var.F"
144	#endif
145
146	c Arguments :
147	c -----------
148
149	c inputs:
150	c -------
151	INTEGER ngrid,nlayer,nq
152	REAL ptimestep
153	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
154	REAL pphi(ngridmx,nlayer)
155	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
156	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
157	REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
158	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
159	LOGICAL firstcall,lastcall
160
161	REAL pday
162	REAL ptime
163	logical tracerdyn
164
165	c outputs:
166	c --------
167	c physical tendencies
168	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
169	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
170	REAL pdpsrf(ngridmx) ! surface pressure tendency
171
172
173	c Local saved variables:
174	c ----------------------
175	c aerosol (dust or ice) extinction optical depth at reference wavelength
176	c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
177	c aerosol optical properties :
178	REAL aerosol(ngridmx,nlayermx,naerkind)
179
180	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
181	INTEGER icount ! counter of calls to physiq during the run.
182	REAL tsurf(ngridmx) ! Surface temperature (K)
183	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
184	REAL co2ice(ngridmx) ! co2 ice surface layer (kg.m-2)
185	REAL albedo(ngridmx,2) ! Surface albedo in each solar band
186	REAL emis(ngridmx) ! Thermal IR surface emissivity
187	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
188	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
189	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
190	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
191	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
192	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
193	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
194
195	c Variables used by the water ice microphysical scheme:
196	REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m)
197	REAL nuice(ngridmx,nlayermx) ! Estimated effective variance
198	! of the size distribution
199	c Albedo of deposited surface ice
200	!!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
201	REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB
202
203	c Variables used by the slope model
204	REAL sl_ls, sl_lct, sl_lat
205	REAL sl_tau, sl_alb, sl_the, sl_psi
206	REAL sl_fl0, sl_flu
207	REAL sl_ra, sl_di0
208	REAL sky
209
210	SAVE day_ini, icount
211	SAVE aerosol, tsurf,tsoil
212	SAVE co2ice,albedo,emis, q2
213	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
214
215	REAL stephan
216	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
217	SAVE stephan
218
219	c Local variables :
220	c -----------------
221
222	REAL CBRT
223	EXTERNAL CBRT
224
225	CHARACTER*80 fichier
226	INTEGER l,ig,ierr,igout,iq,i, tapphys
227
228	REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)
229	REAL fluxsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)
230	REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)
231	REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)
232	REAL tauref(ngridmx) ! Reference column optical depth at 700 Pa
233	! (used if active=F)
234	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
235	REAL zls ! solar longitude (rad)
236	REAL zday ! date (time since Ls=0, in martian days)
237	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
238	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
239	REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
240
241	c Tendancies due to various processes:
242	REAL dqsurf(ngridmx,nqmx)
243	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
244	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
245	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear area
246	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear area
247	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
248	REAL zdtnlte(ngridmx,nlayermx) ! (K/s)
249	REAL zdtsurf(ngridmx) ! (K/s)
250	REAL zdtcloud(ngridmx,nlayermx)
251	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
252	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
253	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
254	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
255	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
256	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
257	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
258	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
259
260	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
261	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
262	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
263	REAL zdqadj(ngridmx,nlayermx,nqmx)
264	REAL zdqc(ngridmx,nlayermx,nqmx)
265	REAL zdqcloud(ngridmx,nlayermx,nqmx)
266	REAL zdqscloud(ngridmx,nqmx)
267	REAL zdqchim(ngridmx,nlayermx,nqmx)
268	REAL zdqschim(ngridmx,nqmx)
269
270	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
271	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
272	REAL zdumolvis(ngridmx,nlayermx)
273	REAL zdvmolvis(ngridmx,nlayermx)
274	real zdqmoldiff(ngridmx,nlayermx,nqmx)
275
276	c Local variable for local intermediate calcul:
277	REAL zflubid(ngridmx)
278	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
279	REAL zdum1(ngridmx,nlayermx)
280	REAL zdum2(ngridmx,nlayermx)
281	REAL ztim1,ztim2,ztim3, z1,z2
282	REAL ztime_fin
283	REAL zdh(ngridmx,nlayermx)
284	INTEGER length
285	PARAMETER (length=100)
286
287	c local variables only used for diagnostic (output in file "diagfi" or "stats")
288	c -----------------------------------------------------------------------------
289	REAL ps(ngridmx), zt(ngridmx,nlayermx)
290	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
291	REAL zq(ngridmx,nlayermx,nqmx)
292	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
293	character*2 str2
294	character*5 str5
295	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
296	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
297	! (particules kg-1)
298	SAVE ccn !! in case iradia != 1
299	real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
300	real qtot1,qtot2 ! total aerosol mass
301	integer igmin, lmin
302	logical tdiag
303
304	real co2col(ngridmx) ! CO2 column
305	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
306	REAL zstress(ngrid), cd
307	real hco2(nqmx),tmean, zlocal(nlayermx)
308	real rho(ngridmx,nlayermx) ! density
309	real vmr(ngridmx,nlayermx) ! volume mixing ratio
310	REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2)
311	REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
312	REAL rave(ngridmx) ! Mean water ice effective radius (m)
313	REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
314	REAL tauTES(ngridmx) ! column optical depth at 825 cm-1
315	REAL Qabsice ! Water ice absorption coefficient
316
317
318	REAL time_phys
319
320	c Variables from thermal
321
322	REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
323	REAL, SAVE :: wmax_th(ngridmx)
324	REAL, SAVE :: hfmax_th(ngridmx)
325	REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
326	REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
327	INTEGER lmax_th(ngridmx)
328	REAL dtke_th(ngridmx,nlayermx+1)
329	REAL zcdv(ngridmx)
330	c=======================================================================
331
332	c 1. Initialisation:
333	c -----------------
334
335	c 1.1 Initialisation only at first call
336	c ---------------------------------------
337	IF (firstcall) THEN
338
339	c variables set to 0
340	c ~~~~~~~~~~~~~~~~~~
341	call zerophys(ngridnlayernaerkind,aerosol)
342	call zerophys(ngrid*nlayer,dtrad)
343	call zerophys(ngrid,fluxrad)
344
345	c read startfi
346	c ~~~~~~~~~~~~
347	#ifndef MESOSCALE
348	! Read netcdf initial physical parameters.
349	CALL phyetat0 ("startfi.nc",0,0,
350	& nsoilmx,nq,
351	& day_ini,time_phys,
352	& tsurf,tsoil,emis,q2,qsurf,co2ice)
353	#else
354	#include "meso_inc/meso_inc_ini.F"
355	#endif
356
357	if (pday.ne.day_ini) then
358	write(,) "PHYSIQ: ERROR: bad synchronization between ",
359	& "physics and dynamics"
360	write(,) "dynamics day: ",pday
361	write(,) "physics day: ",day_ini
362	stop
363	endif
364
365	write (,) 'In physiq day_ini =', day_ini
366
367	c Initialize albedo and orbital calculation
368	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
369	CALL surfini(ngrid,co2ice,qsurf,albedo)
370	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
371
372	c initialize soil
373	c ~~~~~~~~~~~~~~~
374	IF (callsoil) THEN
375	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
376	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
377	ELSE
378	PRINT*,
379	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
380	DO ig=1,ngrid
381	capcal(ig)=1.e5
382	fluxgrd(ig)=0.
383	ENDDO
384	ENDIF
385	icount=1
386
387
388	c initialize tracers
389	c ~~~~~~~~~~~~~~~~~~
390	tracerdyn=tracer
391	IF (tracer) THEN
392	CALL initracer(qsurf,co2ice)
393	ENDIF ! end tracer
394
395	#ifdef MESOSCALE
396	#include "meso_inc/meso_inc_caps.F"
397	#endif
398
399	#ifndef MESOSCALE
400	c Initialize thermospheric parameters
401	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
402
403	if (callthermos) call param_read
404	#endif
405	c Initialize R and Cp as constant
406
407	if (.not.callthermos .and. .not.photochem) then
408	do l=1,nlayermx
409	do ig=1,ngridmx
410	rnew(ig,l)=r
411	cpnew(ig,l)=cpp
412	mmean(ig,l)=mugaz
413	enddo
414	enddo
415	endif
416
417	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
418	write(,)"physiq: water_param Surface ice alb:",alb_surfice
419	ENDIF
420
421	ENDIF ! (end of "if firstcall")
422
423	c ---------------------------------------------------
424	c 1.2 Initializations done at every physical timestep:
425	c ---------------------------------------------------
426	c
427	IF (ngrid.NE.ngridmx) THEN
428	PRINT*,'STOP in PHYSIQ'
429	PRINT*,'Probleme de dimensions :'
430	PRINT*,'ngrid = ',ngrid
431	PRINT*,'ngridmx = ',ngridmx
432	STOP
433	ENDIF
434
435	c Initialize various variables
436	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
437	call zerophys(ngrid*nlayer, pdv)
438	call zerophys(ngrid*nlayer, pdu)
439	call zerophys(ngrid*nlayer, pdt)
440	call zerophys(ngridnlayernq, pdq)
441	call zerophys(ngrid, pdpsrf)
442	call zerophys(ngrid, zflubid)
443	call zerophys(ngrid, zdtsurf)
444	call zerophys(ngrid*nq, dqsurf)
445	igout=ngrid/2+1
446
447
448	zday=pday+ptime ! compute time, in sols (and fraction thereof)
449
450	c Compute Solar Longitude (Ls) :
451	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
452	if (season) then
453	call solarlong(zday,zls)
454	else
455	call solarlong(float(day_ini),zls)
456	end if
457
458	c Compute geopotential at interlayers
459	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
460	c ponderation des altitudes au niveau des couches en dp/p
461
462	DO l=1,nlayer
463	DO ig=1,ngrid
464	zzlay(ig,l)=pphi(ig,l)/g
465	ENDDO
466	ENDDO
467	DO ig=1,ngrid
468	zzlev(ig,1)=0.
469	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
470	ENDDO
471	DO l=2,nlayer
472	DO ig=1,ngrid
473	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
474	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
475	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
476	ENDDO
477	ENDDO
478
479
480	! Potential temperature calculation not the same in physiq and dynamic
481
482	c Compute potential temperature
483	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
484	DO l=1,nlayer
485	DO ig=1,ngrid
486	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
487	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
488	ENDDO
489	ENDDO
490
491	#ifndef MESOSCALE
492	c-----------------------------------------------------------------------
493	c 1.2.5 Compute mean mass, cp, and R
494	c --------------------------------
495
496	if(photochem.or.callthermos) then
497	call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
498	endif
499	#endif
500	c-----------------------------------------------------------------------
501	c 2. Compute radiative tendencies :
502	c------------------------------------
503
504
505	IF (callrad) THEN
506	IF( MOD(icount-1,iradia).EQ.0) THEN
507
508	c Local Solar zenith angle
509	c ~~~~~~~~~~~~~~~~~~~~~~~~
510	CALL orbite(zls,dist_sol,declin)
511
512	IF(diurnal) THEN
513	ztim1=SIN(declin)
514	ztim2=COS(declin)COS(2.pi*(zday-.5))
515	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
516
517	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
518	s ztim1,ztim2,ztim3, mu0,fract)
519
520	ELSE
521	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
522	ENDIF
523
524	c NLTE cooling from CO2 emission
525	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
526
527	IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
528
529	c Find number of layers for LTE radiation calculations
530	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
531	& CALL nlthermeq(ngrid,nlayer,pplev,pplay)
532
533	c Note: Dustopacity.F has been transferred to callradite.F
534
535	c Call main radiative transfer scheme
536	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
537	c Transfer through CO2 (except NIR CO2 absorption)
538	c and aerosols (dust and water ice)
539
540	c Radiative transfer
541	c ------------------
542
543	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
544	$ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
545	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
546	& tauref,tau,aerosol,ccn,rdust,rice,nuice)
547
548	c Outputs for basic check (middle of domain)
549	c ------------------------------------------
550	print, 'check lat lon', lati(igout)180/pi,
551	. long(igout)*180/pi
552	print, 'Ls =',zls180./pi
553	print*, 'tauref(700 Pa) =',tauref(igout)
554	print, 'tau(700 Pa) =',tau(igout,1)700./pplev(igout,1)
555
556	c ---------------------------------------------------------
557	c Call slope parameterization for direct and scattered flux
558	c ---------------------------------------------------------
559	IF(callslope) THEN
560	print *, 'Slope scheme is on and computing...'
561	DO ig=1,ngrid
562	sl_the = theta_sl(ig)
563	IF (sl_the .ne. 0.) THEN
564	ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
565	DO l=1,2
566	sl_lct = ptime24. + 180.long(ig)/pi/15.
567	sl_ra = pi*(1.0-sl_lct/12.)
568	sl_lat = 180.*lati(ig)/pi
569	sl_tau = tau(ig,1)
570	sl_alb = albedo(ig,l)
571	sl_psi = psi_sl(ig)
572	sl_fl0 = fluxsurf_sw(ig,l)
573	sl_di0 = 0.
574	if (mu0(ig) .gt. 0.) then
575	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
576	sl_di0 = sl_di0*1370./dist_sol/dist_sol
577	sl_di0 = sl_di0/ztim1
578	sl_di0 = fluxsurf_sw(ig,l)*sl_di0
579	endif
580	! you never know (roundup concern...)
581	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
582	!!!!!!!!!!!!!!!!!!!!!!!!!!
583	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
584	& sl_tau, sl_alb, sl_the, sl_psi,
585	& sl_di0, sl_fl0, sl_flu )
586	!!!!!!!!!!!!!!!!!!!!!!!!!!
587	fluxsurf_sw(ig,l) = sl_flu
588	ENDDO
589	!!! compute correction on IR flux as well
590	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
591	fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
592	ENDIF
593	ENDDO
594	ENDIF
595
596	c CO2 near infrared absorption
597	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
598	call zerophys(ngrid*nlayer,zdtnirco2)
599	if (callnirco2) then
600	call nirco2abs (ngrid,nlayer,pplay,dist_sol,
601	. mu0,fract,declin, zdtnirco2)
602	endif
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)*(1.-albedo(ig,1))
609	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
610	ENDDO
611
612
613	c Net atmospheric radiative heating rate (K.s-1)
614	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
615	IF(callnlte) THEN
616	CALL blendrad(ngrid, nlayer, pplay,
617	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
618	ELSE
619	DO l=1,nlayer
620	DO ig=1,ngrid
621	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
622	& +zdtnirco2(ig,l)
623	ENDDO
624	ENDDO
625	ENDIF
626
627	ENDIF ! of if(mod(icount-1,iradia).eq.0)
628
629	c Transformation of the radiative tendencies:
630	c -------------------------------------------
631
632	c Net radiative surface flux (W.m-2)
633	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
634	c
635	DO ig=1,ngrid
636	zplanck(ig)=tsurf(ig)*tsurf(ig)
637	zplanck(ig)=emis(ig)*
638	$ stephanzplanck(ig)zplanck(ig)
639	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
640	IF(callslope) THEN
641	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
642	fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
643	ENDIF
644	ENDDO
645
646	DO l=1,nlayer
647	DO ig=1,ngrid
648	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
649	ENDDO
650	ENDDO
651
652	ENDIF ! of IF (callrad)
653
654	c-----------------------------------------------------------------------
655	c 3. Gravity wave and subgrid scale topography drag :
656	c -------------------------------------------------
657
658
659	IF(calllott)THEN
660
661	CALL calldrag_noro(ngrid,nlayer,ptimestep,
662	& pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
663
664	DO l=1,nlayer
665	DO ig=1,ngrid
666	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
667	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
668	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
669	ENDDO
670	ENDDO
671	ENDIF
672
673	c-----------------------------------------------------------------------
674	c 4. Vertical diffusion (turbulent mixing):
675	c -----------------------------------------
676
677	IF (calldifv) THEN
678
679	DO ig=1,ngrid
680	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
681	ENDDO
682
683	CALL zerophys(ngrid*nlayer,zdum1)
684	CALL zerophys(ngrid*nlayer,zdum2)
685	do l=1,nlayer
686	do ig=1,ngrid
687	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
688	enddo
689	enddo
690
691	c Calling vdif (Martian version WITH CO2 condensation)
692	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
693	$ ptimestep,capcal,lwrite,
694	$ pplay,pplev,zzlay,zzlev,z0,
695	$ pu,pv,zh,pq,tsurf,emis,qsurf,
696	$ zdum1,zdum2,zdh,pdq,zflubid,
697	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
698	& zdqdif,zdqsdif,wmax_th,zcdv)
699
700	#ifdef MESOSCALE
701	#include "meso_inc/meso_inc_les.F"
702	#endif
703	DO l=1,nlayer
704	DO ig=1,ngrid
705	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
706	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
707	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
708
709	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
710
711	ENDDO
712	ENDDO
713
714	DO ig=1,ngrid
715	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
716	ENDDO
717
718	if (tracer) then
719	DO iq=1, nq
720	DO l=1,nlayer
721	DO ig=1,ngrid
722	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
723	ENDDO
724	ENDDO
725	ENDDO
726	DO iq=1, nq
727	DO ig=1,ngrid
728	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
729	ENDDO
730	ENDDO
731	end if ! of if (tracer)
732
733	ELSE
734	DO ig=1,ngrid
735	zdtsurf(ig)=zdtsurf(ig)+
736	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
737	ENDDO
738	#ifdef MESOSCALE
739	IF (flag_LES) THEN
740	write(,) 'LES mode !'
741	write(,) 'Please set calldifv to T in callphys.def'
742	STOP
743	ENDIF
744	#endif
745	ENDIF ! of IF (calldifv)
746
747	c-----------------------------------------------------------------------
748	c TEST. Thermals :
749	c HIGHLY EXPERIMENTAL, BEWARE !!
750	c -----------------------------
751
752	if(calltherm) then
753
754	call calltherm_interface(firstcall,
755	$ long,lati,zzlev,zzlay,
756	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
757	$ pplay,pplev,pphi,zpopsk,
758	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
759	$ dtke_th,hfmax_th,wmax_th)
760
761	DO l=1,nlayer
762	DO ig=1,ngrid
763	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
764	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
765	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
766	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
767	ENDDO
768	ENDDO
769
770	DO ig=1,ngrid
771	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
772	ENDDO
773
774	if (tracer) then
775	DO iq=1,nq
776	DO l=1,nlayer
777	DO ig=1,ngrid
778	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
779	ENDDO
780	ENDDO
781	ENDDO
782	endif
783
784	else !of if calltherm
785	lmax_th(:)=0
786	end if
787
788	c-----------------------------------------------------------------------
789	c 5. Dry convective adjustment:
790	c -----------------------------
791
792	IF(calladj) THEN
793
794	DO l=1,nlayer
795	DO ig=1,ngrid
796	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
797	ENDDO
798	ENDDO
799	CALL zerophys(ngrid*nlayer,zduadj)
800	CALL zerophys(ngrid*nlayer,zdvadj)
801	CALL zerophys(ngrid*nlayer,zdhadj)
802
803	CALL convadj(ngrid,nlayer,nq,ptimestep,
804	$ pplay,pplev,zpopsk,lmax_th,
805	$ pu,pv,zh,pq,
806	$ pdu,pdv,zdh,pdq,
807	$ zduadj,zdvadj,zdhadj,
808	$ zdqadj)
809
810
811	DO l=1,nlayer
812	DO ig=1,ngrid
813	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
814	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
815	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
816
817	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
818	ENDDO
819	ENDDO
820
821	if(tracer) then
822	DO iq=1, nq
823	DO l=1,nlayer
824	DO ig=1,ngrid
825	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
826	ENDDO
827	ENDDO
828	ENDDO
829	end if
830	ENDIF ! of IF(calladj)
831
832	c-----------------------------------------------------------------------
833	c 6. Carbon dioxide condensation-sublimation:
834	c -------------------------------------------
835
836	IF (callcond) THEN
837	CALL newcondens(ngrid,nlayer,nq,ptimestep,
838	$ capcal,pplay,pplev,tsurf,pt,
839	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
840	$ co2ice,albedo,emis,
841	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
842	$ fluxsurf_sw,zls)
843
844	DO l=1,nlayer
845	DO ig=1,ngrid
846	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
847	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
848	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
849	ENDDO
850	ENDDO
851	DO ig=1,ngrid
852	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
853	ENDDO
854
855	IF (tracer) THEN
856	DO iq=1, nq
857	DO l=1,nlayer
858	DO ig=1,ngrid
859	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
860	ENDDO
861	ENDDO
862	ENDDO
863	ENDIF ! of IF (tracer)
864
865	ENDIF ! of IF (callcond)
866
867	c-----------------------------------------------------------------------
868	c 7. Specific parameterizations for tracers
869	c: -----------------------------------------
870
871	if (tracer) then
872
873	c 7a. Water and ice
874	c ---------------
875
876	c ---------------------------------------
877	c Water ice condensation in the atmosphere
878	c ----------------------------------------
879	IF (water) THEN
880
881	call watercloud(ngrid,nlayer,ptimestep,
882	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
883	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
884	& nq,naerkind,tau,
885	& ccn,rdust,rice,nuice)
886	if (activice) then
887	c Temperature variation due to latent heat release
888	DO l=1,nlayer
889	DO ig=1,ngrid
890	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
891	ENDDO
892	ENDDO
893	endif
894
895	! increment water vapour and ice atmospheric tracers tendencies
896	IF (water) THEN
897	DO l=1,nlayer
898	DO ig=1,ngrid
899	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
900	& zdqcloud(ig,l,igcm_h2o_vap)
901	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
902	& zdqcloud(ig,l,igcm_h2o_ice)
903	ENDDO
904	ENDDO
905	ENDIF ! of IF (water) THEN
906	! Increment water ice surface tracer tendency
907	DO ig=1,ngrid
908	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
909	& zdqscloud(ig,igcm_h2o_ice)
910	ENDDO
911
912	END IF ! of IF (water)
913
914	c 7b. Chemical species
915	c ------------------
916
917	#ifndef MESOSCALE
918	c --------------
919	c photochemistry :
920	c --------------
921	IF (photochem .or. thermochem) then
922	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
923	& zzlay,zday,pq,pdq,rice,
924	& zdqchim,zdqschim,zdqcloud,zdqscloud)
925	!NB: Photochemistry includes condensation of H2O2
926
927	! increment values of tracers:
928	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
929	! tracers is zero anyways
930	DO l=1,nlayer
931	DO ig=1,ngrid
932	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
933	ENDDO
934	ENDDO
935	ENDDO ! of DO iq=1,nq
936	! add condensation tendency for H2O2
937	if (igcm_h2o2.ne.0) then
938	DO l=1,nlayer
939	DO ig=1,ngrid
940	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
941	& +zdqcloud(ig,l,igcm_h2o2)
942	ENDDO
943	ENDDO
944	endif
945
946	! increment surface values of tracers:
947	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
948	! tracers is zero anyways
949	DO ig=1,ngrid
950	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
951	ENDDO
952	ENDDO ! of DO iq=1,nq
953	! add condensation tendency for H2O2
954	if (igcm_h2o2.ne.0) then
955	DO ig=1,ngrid
956	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
957	& +zdqscloud(ig,igcm_h2o2)
958	ENDDO
959	endif
960
961	END IF ! of IF (photochem.or.thermochem)
962	#endif
963
964	c 7c. Aerosol particles
965	c -------------------
966
967	c ----------
968	c Dust devil :
969	c ----------
970	IF(callddevil) then
971	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
972	& zdqdev,zdqsdev)
973
974	if (dustbin.ge.1) then
975	do iq=1,nq
976	DO l=1,nlayer
977	DO ig=1,ngrid
978	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
979	ENDDO
980	ENDDO
981	enddo
982	do iq=1,nq
983	DO ig=1,ngrid
984	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
985	ENDDO
986	enddo
987	endif ! of if (dustbin.ge.1)
988
989	END IF ! of IF (callddevil)
990
991	c -------------
992	c Sedimentation : acts also on water ice
993	c -------------
994	IF (sedimentation) THEN
995	!call zerophys(ngridnlayernq, zdqsed)
996	zdqsed(1:ngrid,1:nlayer,1:nq)=0
997	!call zerophys(ngrid*nq, zdqssed)
998	zdqssed(1:ngrid,1:nq)=0
999
1000	call callsedim(ngrid,nlayer, ptimestep,
1001	& pplev,zzlev, pt, rdust, rice,
1002	& pq, pdq, zdqsed, zdqssed,nq)
1003
1004	DO iq=1, nq
1005	DO l=1,nlayer
1006	DO ig=1,ngrid
1007	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1008	ENDDO
1009	ENDDO
1010	ENDDO
1011	DO iq=1, nq
1012	DO ig=1,ngrid
1013	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1014	ENDDO
1015	ENDDO
1016	END IF ! of IF (sedimentation)
1017
1018	c 7d. Updates
1019	c ---------
1020
1021	DO iq=1, nq
1022	DO ig=1,ngrid
1023
1024	c ---------------------------------
1025	c Updating tracer budget on surface
1026	c ---------------------------------
1027	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1028
1029	ENDDO ! (ig)
1030	ENDDO ! (iq)
1031
1032	endif ! of if (tracer)
1033
1034	#ifndef MESOSCALE
1035	c-----------------------------------------------------------------------
1036	c 8. THERMOSPHERE CALCULATION
1037	c-----------------------------------------------------------------------
1038
1039	if (callthermos) then
1040	call thermosphere(pplev,pplay,dist_sol,
1041	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1042	& pt,pq,pu,pv,pdt,pdq,
1043	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1044
1045	DO l=1,nlayer
1046	DO ig=1,ngrid
1047	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1048	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1049	& +zdteuv(ig,l)
1050	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1051	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1052	DO iq=1, nq
1053	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1054	ENDDO
1055	ENDDO
1056	ENDDO
1057
1058	endif ! of if (callthermos)
1059	#endif
1060
1061	c-----------------------------------------------------------------------
1062	c 9. Surface and sub-surface soil temperature
1063	c-----------------------------------------------------------------------
1064	c
1065	c
1066	c 9.1 Increment Surface temperature:
1067	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1068
1069	DO ig=1,ngrid
1070	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1071	ENDDO
1072
1073	c Prescribe a cold trap at south pole (except at high obliquity !!)
1074	c Temperature at the surface is set there to be the temperature
1075	c corresponding to equilibrium temperature between phases of CO2
1076
1077	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
1078	#ifndef MESOSCALE
1079	if (caps.and.(obliquit.lt.27.)) then
1080	! NB: Updated surface pressure, at grid point 'ngrid', is
1081	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1082	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1083	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1084	endif
1085	#endif
1086	c -------------------------------------------------------------
1087	c Change of surface albedo (set to 0.4) in case of ground frost
1088	c everywhere except on the north permanent cap and in regions
1089	c covered by dry ice.
1090	c ALWAYS PLACE these lines after newcondens !!!
1091	c -------------------------------------------------------------
1092	do ig=1,ngrid
1093	if ((co2ice(ig).eq.0).and.
1094	& (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
1095	albedo(ig,1) = alb_surfice
1096	albedo(ig,2) = alb_surfice
1097	endif
1098	enddo ! of do ig=1,ngrid
1099	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
1100
1101	c
1102	c 9.2 Compute soil temperatures and subsurface heat flux:
1103	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1104	IF (callsoil) THEN
1105	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1106	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
1107	ENDIF
1108
1109	c-----------------------------------------------------------------------
1110	c 10. Write output files
1111	c ----------------------
1112
1113	c -------------------------------
1114	c Dynamical fields incrementation
1115	c -------------------------------
1116	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1117	! temperature, zonal and meridional wind
1118	DO l=1,nlayer
1119	DO ig=1,ngrid
1120	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1121	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1122	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1123	ENDDO
1124	ENDDO
1125
1126	! tracers
1127	DO iq=1, nq
1128	DO l=1,nlayer
1129	DO ig=1,ngrid
1130	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1131	ENDDO
1132	ENDDO
1133	ENDDO
1134
1135	! surface pressure
1136	DO ig=1,ngrid
1137	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1138	ENDDO
1139
1140	! pressure
1141	DO l=1,nlayer
1142	DO ig=1,ngrid
1143	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1144	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1145	ENDDO
1146	ENDDO
1147
1148	! Density
1149	DO l=1,nlayer
1150	DO ig=1,ngrid
1151	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1152	ENDDO
1153	ENDDO
1154
1155	c Compute surface stress : (NB: z0 is a common in surfdat.h)
1156	c DO ig=1,ngrid
1157	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
1158	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1159	c ENDDO
1160
1161	c Sum of fluxes in solar spectral bands (for output only)
1162	DO ig=1,ngrid
1163	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1164	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1165	ENDDO
1166	c ***** TEST ****************************************************
1167	ztim1 = 999
1168	DO l=1,nlayer
1169	DO ig=1,ngrid
1170	if (pt(ig,l).lt.ztim1) then
1171	ztim1 = pt(ig,l)
1172	igmin = ig
1173	lmin = l
1174	end if
1175	ENDDO
1176	ENDDO
1177	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1178	write(,) 'PHYSIQ: stability WARNING :'
1179	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1180	& 'ig l =', igmin, lmin
1181	end if
1182	c *******************************************************************
1183
1184	c ---------------------
1185	c Outputs to the screen
1186	c ---------------------
1187
1188	IF (lwrite) THEN
1189	PRINT*,'Global diagnostics for the physics'
1190	PRINT,'Variables and their increments x and dx/dt dt'
1191	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1192	WRITE(*,'(2f10.5,2f15.5)')
1193	s tsurf(igout),zdtsurf(igout)*ptimestep,
1194	s pplev(igout,1),pdpsrf(igout)*ptimestep
1195	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1196	WRITE(*,'(i4,6f10.5)') (l,
1197	s pu(igout,l),pdu(igout,l)*ptimestep,
1198	s pv(igout,l),pdv(igout,l)*ptimestep,
1199	s pt(igout,l),pdt(igout,l)*ptimestep,
1200	s l=1,nlayer)
1201	ENDIF ! of IF (lwrite)
1202
1203	IF (ngrid.NE.1) THEN
1204
1205	#ifndef MESOSCALE
1206	c -------------------------------------------------------------------
1207	c Writing NetCDF file "RESTARTFI" at the end of the run
1208	c -------------------------------------------------------------------
1209	c Note: 'restartfi' is stored just before dynamics are stored
1210	c in 'restart'. Between now and the writting of 'restart',
1211	c there will have been the itau=itau+1 instruction and
1212	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1213	c thus we store for time=time+dtvr
1214
1215	IF(lastcall) THEN
1216	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1217	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1218	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1219	. ptimestep,pday,
1220	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
1221	. area,albedodat,inertiedat,zmea,zstd,zsig,
1222	. zgam,zthe)
1223	ENDIF
1224	#endif
1225
1226	c -------------------------------------------------------------------
1227	c Calculation of diagnostic variables written in both stats and
1228	c diagfi files
1229	c -------------------------------------------------------------------
1230
1231	if (tracer) then
1232	if (water) then
1233
1234	call zerophys(ngrid,mtot)
1235	call zerophys(ngrid,icetot)
1236	call zerophys(ngrid,rave)
1237	call zerophys(ngrid,tauTES)
1238	do ig=1,ngrid
1239	do l=1,nlayermx
1240	mtot(ig) = mtot(ig) +
1241	& zq(ig,l,igcm_h2o_vap) *
1242	& (pplev(ig,l) - pplev(ig,l+1)) / g
1243	icetot(ig) = icetot(ig) +
1244	& zq(ig,l,igcm_h2o_ice) *
1245	& (pplev(ig,l) - pplev(ig,l+1)) / g
1246	rave(ig) = rave(ig) +
1247	& zq(ig,l,igcm_h2o_ice) *
1248	& (pplev(ig,l) - pplev(ig,l+1)) / g *
1249	& rice(ig,l) * (1.+nuice_ref)
1250	c Computing abs optical depth at 825 cm-1 in each
1251	c layer to simulate NEW TES retrieval
1252	Qabsice = min(
1253	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
1254	& )
1255	opTES(ig,l)= 0.75 * Qabsice *
1256	& zq(ig,l,igcm_h2o_ice) *
1257	& (pplev(ig,l) - pplev(ig,l+1)) / g
1258	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1259	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1260	enddo
1261	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
1262	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1263	enddo
1264
1265	endif ! of if (water)
1266	endif ! of if (tracer)
1267
1268	c -----------------------------------------------------------------
1269	c WSTATS: Saving statistics
1270	c -----------------------------------------------------------------
1271	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1272	c which can later be used to make the statistic files of the run:
1273	c "stats") only possible in 3D runs !
1274
1275	IF (callstats) THEN
1276
1277	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1278	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1279	call wstats(ngrid,"co2ice","CO2 ice cover",
1280	& "kg.m-2",2,co2ice)
1281	call wstats(ngrid,"fluxsurf_lw",
1282	& "Thermal IR radiative flux to surface","W.m-2",2,
1283	& fluxsurf_lw)
1284	call wstats(ngrid,"fluxsurf_sw",
1285	& "Solar radiative flux to surface","W.m-2",2,
1286	& fluxsurf_sw_tot)
1287	call wstats(ngrid,"fluxtop_lw",
1288	& "Thermal IR radiative flux to space","W.m-2",2,
1289	& fluxtop_lw)
1290	call wstats(ngrid,"fluxtop_sw",
1291	& "Solar radiative flux to space","W.m-2",2,
1292	& fluxtop_sw_tot)
1293	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1294	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1295	call wstats(ngrid,"v","Meridional (North-South) wind",
1296	& "m.s-1",3,zv)
1297	call wstats(ngrid,"w","Vertical (down-up) wind",
1298	& "m.s-1",3,pw)
1299	call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
1300	call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
1301	c call wstats(ngrid,"q2",
1302	c & "Boundary layer eddy kinetic energy",
1303	c & "m2.s-2",3,q2)
1304	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1305	c & emis)
1306	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
1307	c & 2,zstress)
1308	c call wstats(ngrid,"sw_htrt","sw heat.rate",
1309	c & "W.m-2",3,zdtsw)
1310	c call wstats(ngrid,"lw_htrt","lw heat.rate",
1311	c & "W.m-2",3,zdtlw)
1312
1313	if (tracer) then
1314	if (water) then
1315	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1316	& *mugaz/mmol(igcm_h2o_vap)
1317	call wstats(ngrid,"vmr_h2ovapor",
1318	& "H2O vapor volume mixing ratio","mol/mol",
1319	& 3,vmr)
1320	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1321	& *mugaz/mmol(igcm_h2o_ice)
1322	call wstats(ngrid,"vmr_h2oice",
1323	& "H2O ice volume mixing ratio","mol/mol",
1324	& 3,vmr)
1325	call wstats(ngrid,"h2o_ice_s",
1326	& "surface h2o_ice","kg/m2",
1327	& 2,qsurf(1,igcm_h2o_ice))
1328
1329	call wstats(ngrid,"mtot",
1330	& "total mass of water vapor","kg/m2",
1331	& 2,mtot)
1332	call wstats(ngrid,"icetot",
1333	& "total mass of water ice","kg/m2",
1334	& 2,icetot)
1335	call wstats(ngrid,"reffice",
1336	& "Mean reff","m",
1337	& 2,rave)
1338	c call wstats(ngrid,"rice",
1339	c & "Ice particle size","m",
1340	c & 3,rice)
1341	c If activice is true, tauTES is computed in aeropacity.F;
1342	if (.not.activice) then
1343	call wstats(ngrid,"tauTESap",
1344	& "tau abs 825 cm-1","",
1345	& 2,tauTES)
1346	endif
1347
1348	endif ! of if (water)
1349
1350	if (thermochem.or.photochem) then
1351	do iq=1,nq
1352	if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
1353	. (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
1354	. (noms(iq).eq."h2").or.
1355	. (noms(iq).eq."o3")) then
1356	do l=1,nlayer
1357	do ig=1,ngrid
1358	vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
1359	end do
1360	end do
1361	call wstats(ngrid,"vmr_"//trim(noms(iq)),
1362	. "Volume mixing ratio","mol/mol",3,vmr)
1363	endif
1364	enddo
1365	endif ! of if (thermochem.or.photochem)
1366
1367	endif ! of if (tracer)
1368
1369	IF(lastcall) THEN
1370	write (,) "Writing stats..."
1371	call mkstats(ierr)
1372	ENDIF
1373
1374	ENDIF !if callstats
1375
1376	c (Store EOF for Mars Climate database software)
1377	IF (calleofdump) THEN
1378	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1379	ENDIF
1380
1381
1382	#ifdef MESOSCALE
1383	!!!
1384	!!! OUTPUT FIELDS
1385	!!!
1386	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
1387	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
1388	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
1389	IF (igcm_h2o_ice .ne. 0) THEN
1390	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
1391	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1392	. * mugaz / mmol(igcm_h2o_ice)
1393	ENDIF
1394	c AUTOMATICALLY GENERATED FROM REGISTRY
1395	#include "fill_save.inc"
1396	#else
1397
1398	c ==========================================================
1399	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
1400	c any variable for diagnostic (output with period
1401	c "ecritphy", set in "run.def")
1402	c ==========================================================
1403	c WRITEDIAGFI can ALSO be called from any other subroutines
1404	c for any variables !!
1405	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1406	c & emis)
1407	! call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
1408	! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
1409	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1410	& tsurf)
1411	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1412	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
1413	& co2ice)
1414	c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
1415	c & "K",2,zt(1,7))
1416	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1417	& fluxsurf_lw)
1418	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1419	& fluxsurf_sw_tot)
1420	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1421	& fluxtop_lw)
1422	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1423	& fluxtop_sw_tot)
1424	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1425	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1426	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1427	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1428	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1429	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1430	! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1431	c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
1432	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
1433	c & zstress)
1434	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
1435	c & 'w.m-2',3,zdtsw)
1436	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
1437	c & 'w.m-2',3,zdtlw)
1438	c CALL WRITEDIAGFI(ngridmx,'tauTESap',
1439	c & 'tau abs 825 cm-1',
1440	c & '',2,tauTES)
1441
1442
1443	c ----------------------------------------------------------
1444	c Outputs of the CO2 cycle
1445	c ----------------------------------------------------------
1446
1447	if (tracer.and.(igcm_co2.ne.0)) then
1448	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
1449	! & "kg/kg",2,zq(1,1,igcm_co2))
1450	! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
1451	! & "kg/kg",3,zq(1,1,igcm_co2))
1452
1453	! Compute co2 column
1454	call zerophys(ngrid,co2col)
1455	do l=1,nlayermx
1456	do ig=1,ngrid
1457	co2col(ig)=co2col(ig)+
1458	& zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
1459	enddo
1460	enddo
1461	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
1462	& co2col)
1463	endif ! of if (tracer.and.(igcm_co2.ne.0))
1464
1465	c ----------------------------------------------------------
1466	c Outputs of the water cycle
1467	c ----------------------------------------------------------
1468	if (tracer) then
1469	if (water) then
1470
1471	CALL WRITEDIAGFI(ngridmx,'mtot',
1472	& 'total mass of water vapor',
1473	& 'kg/m2',2,mtot)
1474	CALL WRITEDIAGFI(ngridmx,'icetot',
1475	& 'total mass of water ice',
1476	& 'kg/m2',2,icetot)
1477	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1478	c & *mugaz/mmol(igcm_h2o_ice)
1479	c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
1480	c & 'mol/mol',3,vmr)
1481	CALL WRITEDIAGFI(ngridmx,'reffice',
1482	& 'Mean reff',
1483	& 'm',2,rave)
1484	c call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
1485	c & 'm',3,rice)
1486	c If activice is true, tauTES is computed in aeropacity.F;
1487	if (.not.activice) then
1488	CALL WRITEDIAGFI(ngridmx,'tauTESap',
1489	& 'tau abs 825 cm-1',
1490	& '',2,tauTES)
1491	endif
1492	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
1493	& 'surface h2o_ice',
1494	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
1495	endif !(water)
1496
1497
1498	if (water.and..not.photochem) then
1499	iq=nq
1500	c write(str2(1:2),'(i2.2)') iq
1501	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
1502	c & 'kg.m-2',2,zdqscloud(1,iq))
1503	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
1504	c & 'kg/kg',3,zdqchim(1,1,iq))
1505	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
1506	c & 'kg/kg',3,zdqdif(1,1,iq))
1507	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
1508	c & 'kg/kg',3,zdqadj(1,1,iq))
1509	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
1510	c & 'kg/kg',3,zdqc(1,1,iq))
1511	endif !(water.and..not.photochem)
1512	endif
1513
1514	c ----------------------------------------------------------
1515	c Outputs of the dust cycle
1516	c ----------------------------------------------------------
1517
1518	c call WRITEDIAGFI(ngridmx,'tauref',
1519	c & 'Dust ref opt depth','NU',2,tauref)
1520
1521	if (tracer.and.(dustbin.ne.0)) then
1522	c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
1523	if (doubleq) then
1524	c call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
1525	c & 'kg.m-2',2,qsurf(1,1))
1526	c call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
1527	c & 'N.m-2',2,qsurf(1,2))
1528	c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
1529	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
1530	c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
1531	c & 'kg.m-2.s-1',2,zdqssed(1,1))
1532	call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
1533	& 'm',3,rdust*ref_r0)
1534	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
1535	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
1536	c call WRITEDIAGFI(ngridmx,'dustN','Dust number',
1537	c & 'part/kg',3,pq(1,1,igcm_dust_number))
1538	else
1539	do iq=1,dustbin
1540	write(str2(1:2),'(i2.2)') iq
1541	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
1542	& 'kg/kg',3,zq(1,1,iq))
1543	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
1544	& 'kg.m-2',2,qsurf(1,iq))
1545	end do
1546	endif ! (doubleq)
1547	c if (submicron) then
1548	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
1549	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
1550	c endif ! (submicron)
1551	end if ! (tracer.and.(dustbin.ne.0))
1552
1553	c ----------------------------------------------------------
1554	c Outputs of thermals
1555	c ----------------------------------------------------------
1556	if (calltherm) then
1557
1558	! call WRITEDIAGFI(ngrid,'dtke',
1559	! & 'tendance tke thermiques','m2/s2',
1560	! & 3,dtke_th)
1561	! call WRITEDIAGFI(ngrid,'d_u_ajs',
1562	! & 'tendance u thermiques','m/s',
1563	! & 3,pdu_th*ptimestep)
1564	! call WRITEDIAGFI(ngrid,'d_v_ajs',
1565	! & 'tendance v thermiques','m/s',
1566	! & 3,pdv_th*ptimestep)
1567	! if (tracer) then
1568	! if (nq .eq. 2) then
1569	! call WRITEDIAGFI(ngrid,'deltaq_th',
1570	! & 'delta q thermiques','kg/kg',
1571	! & 3,ptimestep*pdq_th(:,:,2))
1572	! endif
1573	! endif
1574
1575	lmax_th_out(:)=real(lmax_th(:))
1576
1577	call WRITEDIAGFI(ngridmx,'lmax_th',
1578	& 'hauteur du thermique','K',
1579	& 2,lmax_th_out)
1580	call WRITEDIAGFI(ngridmx,'hfmax_th',
1581	& 'maximum TH heat flux','K.m/s',
1582	& 2,hfmax_th)
1583	call WRITEDIAGFI(ngridmx,'wmax_th',
1584	& 'maximum TH vertical velocity','m/s',
1585	& 2,wmax_th)
1586
1587	endif
1588
1589	c ----------------------------------------------------------
1590	c Output in netcdf file "diagsoil.nc" for subterranean
1591	c variables (output every "ecritphy", as for writediagfi)
1592	c ----------------------------------------------------------
1593
1594	! Write soil temperature
1595	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
1596	! & 3,tsoil)
1597	! Write surface temperature
1598	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
1599	! & 2,tsurf)
1600
1601	c ==========================================================
1602	c END OF WRITEDIAGFI
1603	c ==========================================================
1604	#endif
1605
1606	ELSE ! if(ngrid.eq.1)
1607
1608	print,'Ls =',zls180./pi,
1609	& ' tauref(700 Pa) =',tauref
1610	c ----------------------------------------------------------------------
1611	c Output in grads file "g1d" (ONLY when using testphys1d)
1612	c (output at every X physical timestep)
1613	c ----------------------------------------------------------------------
1614	c
1615	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
1616	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
1617	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
1618
1619	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
1620	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
1621	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
1622	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
1623
1624	! THERMALS STUFF 1D
1625
1626	if(calltherm) then
1627
1628	lmax_th_out(:)=real(lmax_th(:))
1629
1630	call WRITEDIAGFI(ngridmx,'lmax_th',
1631	& 'hauteur du thermique','point',
1632	& 0,lmax_th_out)
1633	call WRITEDIAGFI(ngridmx,'hfmax_th',
1634	& 'maximum TH heat flux','K.m/s',
1635	& 0,hfmax_th)
1636	call WRITEDIAGFI(ngridmx,'wmax_th',
1637	& 'maximum TH vertical velocity','m/s',
1638	& 0,wmax_th)
1639
1640
1641	co2col(:)=0.
1642	if (tracer) then
1643	do l=1,nlayermx
1644	do ig=1,ngrid
1645	co2col(ig)=co2col(ig)+
1646	& zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
1647	enddo
1648	enddo
1649
1650	end if
1651	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
1652	& ,'kg/m-2',0,co2col)
1653	endif
1654	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
1655	& ,'m/s',1,pw)
1656	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
1657	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
1658	& tsurf)
1659
1660	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
1661	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
1662	! or output in diagfi.nc (for testphys1d)
1663	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
1664	call WRITEDIAGFI(ngridmx,'temp','Temperature',
1665	& 'K',1,zt)
1666
1667	if(tracer) then
1668	c CALL writeg1d(ngrid,1,tau,'tau','SI')
1669	do iq=1,nq
1670	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
1671	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
1672	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
1673	end do
1674	end if
1675
1676	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
1677
1678	do l=2,nlayer-1
1679	tmean=zt(1,l)
1680	if(zt(1,l).ne.zt(1,l-1))
1681	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
1682	zlocal(l)= zlocal(l-1)
1683	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
1684	enddo
1685	zlocal(nlayer)= zlocal(nlayer-1)-
1686	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
1687	& rnew(1,nlayer)*tmean/g
1688
1689	END IF ! if(ngrid.ne.1)
1690
1691	icount=icount+1
1692	RETURN
1693	END

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

Download in other formats:

Original Format