Context Navigation

physiq.F @ 1242

Last change on this file since 1242 was 234, checked in by aslmd, 13 years ago

LMDZ.MARS
+ MESOSCALE

An important change to merge physiq.F and inifis.F for GCM and mesoscale
This is mostly transparent to GCM users and developers (use of MESOSCALE precompiler flags)
Makes it easy (and mostly automatic!) for changes in GCM physics to be impacted in mesoscale physics
A few minor changes have followed in the GCM (slope scheme, one-point diagnostic).

Compilation + run is OK on both sides (GCM and mesoscale).
On the mesoscale side, call_meso_physiq?.inc and call_meso_inifis?.inc have been changed accordingly.

Here is the excerpt from README file:

19/07/2011 == AS

Finished converging meso_physiq.F and meso_inifis.F towards physiq.F and inifis.F --> see previous point 15/07/2011 --> meso_ routines no longer exist (everything is in meso_inc and transparent to GCM users) --> GCM routines include mesoscale parts within MESOSCALE precompiler commands --> MESOSCALE parts are as hidden as possible not to mess up with GCM users/developers
Cleaned inelegant or useless #ifdef [or] #ifndef MESOSCALE in physiq and inifis so that a minimum amount of such precompiler commands is now reached [mainly related to I/O]
Added the SF08 slope insolation model in the general physics parameterizations. Added a callslope keyword in inifis.F and callkeys.h --> This keyword is False by default / True if you use -DMESOSCALE
Removed the obsolete call to Viking Lander 1 diagnostic Replaced it with a diagnostic for opacity at the domain center [valid for GCM and mesoscale]

File size: 61.6 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 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 dummycol(ngridmx)
330
331	c=======================================================================
332
333	c 1. Initialisation:
334	c -----------------
335
336	c 1.1 Initialisation only at first call
337	c ---------------------------------------
338	IF (firstcall) THEN
339
340	c variables set to 0
341	c ~~~~~~~~~~~~~~~~~~
342	call zerophys(ngridnlayernaerkind,aerosol)
343	call zerophys(ngrid*nlayer,dtrad)
344	call zerophys(ngrid,fluxrad)
345
346	c read startfi
347	c ~~~~~~~~~~~~
348	#ifndef MESOSCALE
349	! Read netcdf initial physical parameters.
350	CALL phyetat0 ("startfi.nc",0,0,
351	& nsoilmx,nq,
352	& day_ini,time_phys,
353	& tsurf,tsoil,emis,q2,qsurf,co2ice)
354	#else
355	#include "meso_inc/meso_inc_ini.F"
356	#endif
357
358	if (pday.ne.day_ini) then
359	write(,) "PHYSIQ: ERROR: bad synchronization between ",
360	& "physics and dynamics"
361	write(,) "dynamics day: ",pday
362	write(,) "physics day: ",day_ini
363	stop
364	endif
365
366	write (,) 'In physiq day_ini =', day_ini
367
368	c Initialize albedo and orbital calculation
369	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
370	CALL surfini(ngrid,co2ice,qsurf,albedo)
371	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
372
373	c initialize soil
374	c ~~~~~~~~~~~~~~~
375	IF (callsoil) THEN
376	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
377	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
378	ELSE
379	PRINT*,
380	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
381	DO ig=1,ngrid
382	capcal(ig)=1.e5
383	fluxgrd(ig)=0.
384	ENDDO
385	ENDIF
386	icount=1
387
388
389	c initialize tracers
390	c ~~~~~~~~~~~~~~~~~~
391	tracerdyn=tracer
392	IF (tracer) THEN
393	CALL initracer(qsurf,co2ice)
394	ENDIF ! end tracer
395
396	#ifdef MESOSCALE
397	#include "meso_inc/meso_inc_caps.F"
398	#endif
399
400	#ifndef MESOSCALE
401	c Initialize thermospheric parameters
402	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
403
404	if (callthermos) call param_read
405	#endif
406	c Initialize R and Cp as constant
407
408	if (.not.callthermos .and. .not.photochem) then
409	do l=1,nlayermx
410	do ig=1,ngridmx
411	rnew(ig,l)=r
412	cpnew(ig,l)=cpp
413	mmean(ig,l)=mugaz
414	enddo
415	enddo
416	endif
417
418	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
419	write(,)"physiq: water_param Surface ice alb:",alb_surfice
420	ENDIF
421
422	ENDIF ! (end of "if firstcall")
423
424	c ---------------------------------------------------
425	c 1.2 Initializations done at every physical timestep:
426	c ---------------------------------------------------
427	c
428	IF (ngrid.NE.ngridmx) THEN
429	PRINT*,'STOP in PHYSIQ'
430	PRINT*,'Probleme de dimensions :'
431	PRINT*,'ngrid = ',ngrid
432	PRINT*,'ngridmx = ',ngridmx
433	STOP
434	ENDIF
435
436	c Initialize various variables
437	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
438	call zerophys(ngrid*nlayer, pdv)
439	call zerophys(ngrid*nlayer, pdu)
440	call zerophys(ngrid*nlayer, pdt)
441	call zerophys(ngridnlayernq, pdq)
442	call zerophys(ngrid, pdpsrf)
443	call zerophys(ngrid, zflubid)
444	call zerophys(ngrid, zdtsurf)
445	call zerophys(ngrid*nq, dqsurf)
446	igout=ngrid/2+1
447
448
449	zday=pday+ptime ! compute time, in sols (and fraction thereof)
450
451	c Compute Solar Longitude (Ls) :
452	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
453	if (season) then
454	call solarlong(zday,zls)
455	else
456	call solarlong(float(day_ini),zls)
457	end if
458
459	c Compute geopotential at interlayers
460	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
461	c ponderation des altitudes au niveau des couches en dp/p
462
463	DO l=1,nlayer
464	DO ig=1,ngrid
465	zzlay(ig,l)=pphi(ig,l)/g
466	ENDDO
467	ENDDO
468	DO ig=1,ngrid
469	zzlev(ig,1)=0.
470	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
471	ENDDO
472	DO l=2,nlayer
473	DO ig=1,ngrid
474	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
475	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
476	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
477	ENDDO
478	ENDDO
479
480
481	! Potential temperature calculation not the same in physiq and dynamic
482
483	c Compute potential temperature
484	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
485	DO l=1,nlayer
486	DO ig=1,ngrid
487	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
488	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
489	ENDDO
490	ENDDO
491
492	#ifndef MESOSCALE
493	c-----------------------------------------------------------------------
494	c 1.2.5 Compute mean mass, cp, and R
495	c --------------------------------
496
497	if(photochem.or.callthermos) then
498	call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
499	endif
500	#endif
501	c-----------------------------------------------------------------------
502	c 2. Compute radiative tendencies :
503	c------------------------------------
504
505
506	IF (callrad) THEN
507	IF( MOD(icount-1,iradia).EQ.0) THEN
508
509	c Local Solar zenith angle
510	c ~~~~~~~~~~~~~~~~~~~~~~~~
511	CALL orbite(zls,dist_sol,declin)
512
513	IF(diurnal) THEN
514	ztim1=SIN(declin)
515	ztim2=COS(declin)COS(2.pi*(zday-.5))
516	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
517
518	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
519	s ztim1,ztim2,ztim3, mu0,fract)
520
521	ELSE
522	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
523	ENDIF
524
525	c NLTE cooling from CO2 emission
526	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
527
528	IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
529
530	c Find number of layers for LTE radiation calculations
531	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
532	& CALL nlthermeq(ngrid,nlayer,pplev,pplay)
533
534	c Note: Dustopacity.F has been transferred to callradite.F
535
536	c Call main radiative transfer scheme
537	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
538	c Transfer through CO2 (except NIR CO2 absorption)
539	c and aerosols (dust and water ice)
540
541	c Radiative transfer
542	c ------------------
543
544	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
545	$ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
546	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
547	& tauref,tau,aerosol,ccn,rdust,rice,nuice)
548
549	c Outputs for basic check (middle of domain)
550	c ------------------------------------------
551	print, 'check lat lon', lati(igout)180/pi,
552	. long(igout)*180/pi
553	print, 'Ls =',zls180./pi
554	print*, 'tauref(700 Pa) =',tauref(igout)
555	print, 'tau(700 Pa) =',tau(igout,1)700./pplev(igout,1)
556
557	c ---------------------------------------------------------
558	c Call slope parameterization for direct and scattered flux
559	c ---------------------------------------------------------
560	IF(callslope) THEN
561	print *, 'Slope scheme is on and computing...'
562	DO ig=1,ngrid
563	sl_the = theta_sl(ig)
564	IF (sl_the .ne. 0.) THEN
565	ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
566	DO l=1,2
567	sl_lct = ptime24. + 180.long(ig)/pi/15.
568	sl_ra = pi*(1.0-sl_lct/12.)
569	sl_lat = 180.*lati(ig)/pi
570	sl_tau = tau(ig,1)
571	sl_alb = albedo(ig,l)
572	sl_psi = psi_sl(ig)
573	sl_fl0 = fluxsurf_sw(ig,l)
574	sl_di0 = 0.
575	if (mu0(ig) .gt. 0.) then
576	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
577	sl_di0 = sl_di0*1370./dist_sol/dist_sol
578	sl_di0 = sl_di0/ztim1
579	sl_di0 = fluxsurf_sw(ig,l)*sl_di0
580	endif
581	! you never know (roundup concern...)
582	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
583	!!!!!!!!!!!!!!!!!!!!!!!!!!
584	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
585	& sl_tau, sl_alb, sl_the, sl_psi,
586	& sl_di0, sl_fl0, sl_flu )
587	!!!!!!!!!!!!!!!!!!!!!!!!!!
588	fluxsurf_sw(ig,l) = sl_flu
589	ENDDO
590	!!! compute correction on IR flux as well
591	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
592	fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
593	ENDIF
594	ENDDO
595	ENDIF
596
597	c CO2 near infrared absorption
598	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
599	call zerophys(ngrid*nlayer,zdtnirco2)
600	if (callnirco2) then
601	call nirco2abs (ngrid,nlayer,pplay,dist_sol,
602	. mu0,fract,declin, zdtnirco2)
603	endif
604
605	c Radiative flux from the sky absorbed by the surface (W.m-2)
606	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
607	DO ig=1,ngrid
608	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
609	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
610	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
611	ENDDO
612
613
614	c Net atmospheric radiative heating rate (K.s-1)
615	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
616	IF(callnlte) THEN
617	CALL blendrad(ngrid, nlayer, pplay,
618	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
619	ELSE
620	DO l=1,nlayer
621	DO ig=1,ngrid
622	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
623	& +zdtnirco2(ig,l)
624	ENDDO
625	ENDDO
626	ENDIF
627
628	ENDIF ! of if(mod(icount-1,iradia).eq.0)
629
630	c Transformation of the radiative tendencies:
631	c -------------------------------------------
632
633	c Net radiative surface flux (W.m-2)
634	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
635	c
636	DO ig=1,ngrid
637	zplanck(ig)=tsurf(ig)*tsurf(ig)
638	zplanck(ig)=emis(ig)*
639	$ stephanzplanck(ig)zplanck(ig)
640	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
641	IF(callslope) THEN
642	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
643	fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
644	ENDIF
645	ENDDO
646
647	DO l=1,nlayer
648	DO ig=1,ngrid
649	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
650	ENDDO
651	ENDDO
652
653	ENDIF ! of IF (callrad)
654
655	c-----------------------------------------------------------------------
656	c 3. Gravity wave and subgrid scale topography drag :
657	c -------------------------------------------------
658
659
660	IF(calllott)THEN
661
662	CALL calldrag_noro(ngrid,nlayer,ptimestep,
663	& pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
664
665	DO l=1,nlayer
666	DO ig=1,ngrid
667	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
668	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
669	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
670	ENDDO
671	ENDDO
672	ENDIF
673
674	c-----------------------------------------------------------------------
675	c 4. Vertical diffusion (turbulent mixing):
676	c -----------------------------------------
677
678	IF (calldifv) THEN
679
680	DO ig=1,ngrid
681	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
682	ENDDO
683
684	CALL zerophys(ngrid*nlayer,zdum1)
685	CALL zerophys(ngrid*nlayer,zdum2)
686	do l=1,nlayer
687	do ig=1,ngrid
688	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
689	enddo
690	enddo
691
692	c Calling vdif (Martian version WITH CO2 condensation)
693	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
694	$ ptimestep,capcal,lwrite,
695	$ pplay,pplev,zzlay,zzlev,z0,
696	$ pu,pv,zh,pq,tsurf,emis,qsurf,
697	$ zdum1,zdum2,zdh,pdq,zflubid,
698	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
699	& zdqdif,zdqsdif)
700
701	#ifdef MESOSCALE
702	#include "meso_inc/meso_inc_les.F"
703	#endif
704	DO l=1,nlayer
705	DO ig=1,ngrid
706	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
707	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
708	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
709
710	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
711
712	ENDDO
713	ENDDO
714
715	DO ig=1,ngrid
716	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
717	ENDDO
718
719	if (tracer) then
720	DO iq=1, nq
721	DO l=1,nlayer
722	DO ig=1,ngrid
723	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
724	ENDDO
725	ENDDO
726	ENDDO
727	DO iq=1, nq
728	DO ig=1,ngrid
729	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
730	ENDDO
731	ENDDO
732	end if ! of if (tracer)
733
734	ELSE
735	DO ig=1,ngrid
736	zdtsurf(ig)=zdtsurf(ig)+
737	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
738	ENDDO
739	#ifdef MESOSCALE
740	IF (flag_LES) THEN
741	write(,) 'LES mode !'
742	write(,) 'Please set calldifv to T in callphys.def'
743	STOP
744	ENDIF
745	#endif
746	ENDIF ! of IF (calldifv)
747
748	c-----------------------------------------------------------------------
749	c TEST. Thermals :
750	c HIGHLY EXPERIMENTAL, BEWARE !!
751	c -----------------------------
752
753	if(calltherm) then
754
755	call calltherm_interface(firstcall,
756	$ long,lati,zzlev,zzlay,
757	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
758	$ pplay,pplev,pphi,zpopsk,
759	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
760	$ dtke_th,hfmax_th,wmax_th)
761
762	DO l=1,nlayer
763	DO ig=1,ngrid
764	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
765	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
766	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
767	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
768	ENDDO
769	ENDDO
770
771	DO ig=1,ngrid
772	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
773	ENDDO
774
775	if (tracer) then
776	DO iq=1,nq
777	DO l=1,nlayer
778	DO ig=1,ngrid
779	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
780	ENDDO
781	ENDDO
782	ENDDO
783	endif
784
785	else !of if calltherm
786	lmax_th(:)=0
787	end if
788
789	c-----------------------------------------------------------------------
790	c 5. Dry convective adjustment:
791	c -----------------------------
792
793	IF(calladj) THEN
794
795	DO l=1,nlayer
796	DO ig=1,ngrid
797	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
798	ENDDO
799	ENDDO
800	CALL zerophys(ngrid*nlayer,zduadj)
801	CALL zerophys(ngrid*nlayer,zdvadj)
802	CALL zerophys(ngrid*nlayer,zdhadj)
803
804	CALL convadj(ngrid,nlayer,nq,ptimestep,
805	$ pplay,pplev,zpopsk,lmax_th,
806	$ pu,pv,zh,pq,
807	$ pdu,pdv,zdh,pdq,
808	$ zduadj,zdvadj,zdhadj,
809	$ zdqadj)
810
811
812	DO l=1,nlayer
813	DO ig=1,ngrid
814	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
815	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
816	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
817
818	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
819	ENDDO
820	ENDDO
821
822	if(tracer) then
823	DO iq=1, nq
824	DO l=1,nlayer
825	DO ig=1,ngrid
826	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
827	ENDDO
828	ENDDO
829	ENDDO
830	end if
831	ENDIF ! of IF(calladj)
832
833	c-----------------------------------------------------------------------
834	c 6. Carbon dioxide condensation-sublimation:
835	c -------------------------------------------
836
837	IF (callcond) THEN
838	CALL newcondens(ngrid,nlayer,nq,ptimestep,
839	$ capcal,pplay,pplev,tsurf,pt,
840	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
841	$ co2ice,albedo,emis,
842	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
843	$ fluxsurf_sw,zls)
844
845	DO l=1,nlayer
846	DO ig=1,ngrid
847	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
848	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
849	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
850	ENDDO
851	ENDDO
852	DO ig=1,ngrid
853	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
854	ENDDO
855
856	IF (tracer) THEN
857	DO iq=1, nq
858	DO l=1,nlayer
859	DO ig=1,ngrid
860	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
861	ENDDO
862	ENDDO
863	ENDDO
864	ENDIF ! of IF (tracer)
865
866	ENDIF ! of IF (callcond)
867
868	c-----------------------------------------------------------------------
869	c 7. Specific parameterizations for tracers
870	c: -----------------------------------------
871
872	if (tracer) then
873
874	c 7a. Water and ice
875	c ---------------
876
877	c ---------------------------------------
878	c Water ice condensation in the atmosphere
879	c ----------------------------------------
880	IF (water) THEN
881
882	call watercloud(ngrid,nlayer,ptimestep,
883	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
884	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
885	& nq,naerkind,tau,
886	& ccn,rdust,rice,nuice)
887	if (activice) then
888	c Temperature variation due to latent heat release
889	DO l=1,nlayer
890	DO ig=1,ngrid
891	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
892	ENDDO
893	ENDDO
894	endif
895
896	! increment water vapour and ice atmospheric tracers tendencies
897	IF (water) THEN
898	DO l=1,nlayer
899	DO ig=1,ngrid
900	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
901	& zdqcloud(ig,l,igcm_h2o_vap)
902	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
903	& zdqcloud(ig,l,igcm_h2o_ice)
904	ENDDO
905	ENDDO
906	ENDIF ! of IF (water) THEN
907	! Increment water ice surface tracer tendency
908	DO ig=1,ngrid
909	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
910	& zdqscloud(ig,igcm_h2o_ice)
911	ENDDO
912
913	END IF ! of IF (water)
914
915	c 7b. Chemical species
916	c ------------------
917
918	#ifndef MESOSCALE
919	c --------------
920	c photochemistry :
921	c --------------
922	IF (photochem .or. thermochem) then
923	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
924	& zzlay,zday,pq,pdq,rice,
925	& zdqchim,zdqschim,zdqcloud,zdqscloud)
926	!NB: Photochemistry includes condensation of H2O2
927
928	! increment values of tracers:
929	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
930	! tracers is zero anyways
931	DO l=1,nlayer
932	DO ig=1,ngrid
933	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
934	ENDDO
935	ENDDO
936	ENDDO ! of DO iq=1,nq
937	! add condensation tendency for H2O2
938	if (igcm_h2o2.ne.0) then
939	DO l=1,nlayer
940	DO ig=1,ngrid
941	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
942	& +zdqcloud(ig,l,igcm_h2o2)
943	ENDDO
944	ENDDO
945	endif
946
947	! increment surface values of tracers:
948	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
949	! tracers is zero anyways
950	DO ig=1,ngrid
951	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
952	ENDDO
953	ENDDO ! of DO iq=1,nq
954	! add condensation tendency for H2O2
955	if (igcm_h2o2.ne.0) then
956	DO ig=1,ngrid
957	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
958	& +zdqscloud(ig,igcm_h2o2)
959	ENDDO
960	endif
961
962	END IF ! of IF (photochem.or.thermochem)
963	#endif
964
965	c 7c. Aerosol particles
966	c -------------------
967
968	c ----------
969	c Dust devil :
970	c ----------
971	IF(callddevil) then
972	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
973	& zdqdev,zdqsdev)
974
975	if (dustbin.ge.1) then
976	do iq=1,nq
977	DO l=1,nlayer
978	DO ig=1,ngrid
979	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
980	ENDDO
981	ENDDO
982	enddo
983	do iq=1,nq
984	DO ig=1,ngrid
985	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
986	ENDDO
987	enddo
988	endif ! of if (dustbin.ge.1)
989
990	END IF ! of IF (callddevil)
991
992	c -------------
993	c Sedimentation : acts also on water ice
994	c -------------
995	IF (sedimentation) THEN
996	!call zerophys(ngridnlayernq, zdqsed)
997	zdqsed(1:ngrid,1:nlayer,1:nq)=0
998	!call zerophys(ngrid*nq, zdqssed)
999	zdqssed(1:ngrid,1:nq)=0
1000
1001	call callsedim(ngrid,nlayer, ptimestep,
1002	& pplev,zzlev, pt, rdust, rice,
1003	& pq, pdq, zdqsed, zdqssed,nq)
1004
1005	DO iq=1, nq
1006	DO l=1,nlayer
1007	DO ig=1,ngrid
1008	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1009	ENDDO
1010	ENDDO
1011	ENDDO
1012	DO iq=1, nq
1013	DO ig=1,ngrid
1014	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1015	ENDDO
1016	ENDDO
1017	END IF ! of IF (sedimentation)
1018
1019	c 7d. Updates
1020	c ---------
1021
1022	DO iq=1, nq
1023	DO ig=1,ngrid
1024
1025	c ---------------------------------
1026	c Updating tracer budget on surface
1027	c ---------------------------------
1028	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1029
1030	ENDDO ! (ig)
1031	ENDDO ! (iq)
1032
1033	endif ! of if (tracer)
1034
1035	#ifndef MESOSCALE
1036	c-----------------------------------------------------------------------
1037	c 8. THERMOSPHERE CALCULATION
1038	c-----------------------------------------------------------------------
1039
1040	if (callthermos) then
1041	call thermosphere(pplev,pplay,dist_sol,
1042	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1043	& pt,pq,pu,pv,pdt,pdq,
1044	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1045
1046	DO l=1,nlayer
1047	DO ig=1,ngrid
1048	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1049	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1050	& +zdteuv(ig,l)
1051	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1052	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1053	DO iq=1, nq
1054	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1055	ENDDO
1056	ENDDO
1057	ENDDO
1058
1059	endif ! of if (callthermos)
1060	#endif
1061
1062	c-----------------------------------------------------------------------
1063	c 9. Surface and sub-surface soil temperature
1064	c-----------------------------------------------------------------------
1065	c
1066	c
1067	c 9.1 Increment Surface temperature:
1068	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1069
1070	DO ig=1,ngrid
1071	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1072	ENDDO
1073
1074	c Prescribe a cold trap at south pole (except at high obliquity !!)
1075	c Temperature at the surface is set there to be the temperature
1076	c corresponding to equilibrium temperature between phases of CO2
1077
1078	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
1079	#ifndef MESOSCALE
1080	if (caps.and.(obliquit.lt.27.)) then
1081	! NB: Updated surface pressure, at grid point 'ngrid', is
1082	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1083	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1084	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1085	endif
1086	#endif
1087	c -------------------------------------------------------------
1088	c Change of surface albedo (set to 0.4) in case of ground frost
1089	c everywhere except on the north permanent cap and in regions
1090	c covered by dry ice.
1091	c ALWAYS PLACE these lines after newcondens !!!
1092	c -------------------------------------------------------------
1093	do ig=1,ngrid
1094	if ((co2ice(ig).eq.0).and.
1095	& (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
1096	albedo(ig,1) = alb_surfice
1097	albedo(ig,2) = alb_surfice
1098	endif
1099	enddo ! of do ig=1,ngrid
1100	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
1101
1102	c
1103	c 9.2 Compute soil temperatures and subsurface heat flux:
1104	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1105	IF (callsoil) THEN
1106	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1107	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
1108	ENDIF
1109
1110	c-----------------------------------------------------------------------
1111	c 10. Write output files
1112	c ----------------------
1113
1114	c -------------------------------
1115	c Dynamical fields incrementation
1116	c -------------------------------
1117	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1118	! temperature, zonal and meridional wind
1119	DO l=1,nlayer
1120	DO ig=1,ngrid
1121	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1122	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1123	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1124	ENDDO
1125	ENDDO
1126
1127	! tracers
1128	DO iq=1, nq
1129	DO l=1,nlayer
1130	DO ig=1,ngrid
1131	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1132	ENDDO
1133	ENDDO
1134	ENDDO
1135
1136	! surface pressure
1137	DO ig=1,ngrid
1138	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1139	ENDDO
1140
1141	! pressure
1142	DO l=1,nlayer
1143	DO ig=1,ngrid
1144	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1145	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1146	ENDDO
1147	ENDDO
1148
1149	! Density
1150	DO l=1,nlayer
1151	DO ig=1,ngrid
1152	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1153	ENDDO
1154	ENDDO
1155
1156	c Compute surface stress : (NB: z0 is a common in surfdat.h)
1157	c DO ig=1,ngrid
1158	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
1159	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1160	c ENDDO
1161
1162	c Sum of fluxes in solar spectral bands (for output only)
1163	DO ig=1,ngrid
1164	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1165	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1166	ENDDO
1167	c ***** TEST ****************************************************
1168	ztim1 = 999
1169	DO l=1,nlayer
1170	DO ig=1,ngrid
1171	if (pt(ig,l).lt.ztim1) then
1172	ztim1 = pt(ig,l)
1173	igmin = ig
1174	lmin = l
1175	end if
1176	ENDDO
1177	ENDDO
1178	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1179	write(,) 'PHYSIQ: stability WARNING :'
1180	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1181	& 'ig l =', igmin, lmin
1182	end if
1183	c *******************************************************************
1184
1185	c ---------------------
1186	c Outputs to the screen
1187	c ---------------------
1188
1189	IF (lwrite) THEN
1190	PRINT*,'Global diagnostics for the physics'
1191	PRINT,'Variables and their increments x and dx/dt dt'
1192	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1193	WRITE(*,'(2f10.5,2f15.5)')
1194	s tsurf(igout),zdtsurf(igout)*ptimestep,
1195	s pplev(igout,1),pdpsrf(igout)*ptimestep
1196	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1197	WRITE(*,'(i4,6f10.5)') (l,
1198	s pu(igout,l),pdu(igout,l)*ptimestep,
1199	s pv(igout,l),pdv(igout,l)*ptimestep,
1200	s pt(igout,l),pdt(igout,l)*ptimestep,
1201	s l=1,nlayer)
1202	ENDIF ! of IF (lwrite)
1203
1204	IF (ngrid.NE.1) THEN
1205
1206	#ifndef MESOSCALE
1207	c -------------------------------------------------------------------
1208	c Writing NetCDF file "RESTARTFI" at the end of the run
1209	c -------------------------------------------------------------------
1210	c Note: 'restartfi' is stored just before dynamics are stored
1211	c in 'restart'. Between now and the writting of 'restart',
1212	c there will have been the itau=itau+1 instruction and
1213	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1214	c thus we store for time=time+dtvr
1215
1216	IF(lastcall) THEN
1217	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1218	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1219	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1220	. ptimestep,pday,
1221	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
1222	. area,albedodat,inertiedat,zmea,zstd,zsig,
1223	. zgam,zthe)
1224	ENDIF
1225	#endif
1226
1227	c -------------------------------------------------------------------
1228	c Calculation of diagnostic variables written in both stats and
1229	c diagfi files
1230	c -------------------------------------------------------------------
1231
1232	if (tracer) then
1233	if (water) then
1234
1235	call zerophys(ngrid,mtot)
1236	call zerophys(ngrid,icetot)
1237	call zerophys(ngrid,rave)
1238	call zerophys(ngrid,tauTES)
1239	do ig=1,ngrid
1240	do l=1,nlayermx
1241	mtot(ig) = mtot(ig) +
1242	& zq(ig,l,igcm_h2o_vap) *
1243	& (pplev(ig,l) - pplev(ig,l+1)) / g
1244	icetot(ig) = icetot(ig) +
1245	& zq(ig,l,igcm_h2o_ice) *
1246	& (pplev(ig,l) - pplev(ig,l+1)) / g
1247	rave(ig) = rave(ig) +
1248	& zq(ig,l,igcm_h2o_ice) *
1249	& (pplev(ig,l) - pplev(ig,l+1)) / g *
1250	& rice(ig,l) * (1.+nuice_ref)
1251	c Computing abs optical depth at 825 cm-1 in each
1252	c layer to simulate NEW TES retrieval
1253	Qabsice = min(
1254	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
1255	& )
1256	opTES(ig,l)= 0.75 * Qabsice *
1257	& zq(ig,l,igcm_h2o_ice) *
1258	& (pplev(ig,l) - pplev(ig,l+1)) / g
1259	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1260	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1261	enddo
1262	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
1263	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1264	enddo
1265
1266	endif ! of if (water)
1267	endif ! of if (tracer)
1268
1269	c -----------------------------------------------------------------
1270	c WSTATS: Saving statistics
1271	c -----------------------------------------------------------------
1272	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1273	c which can later be used to make the statistic files of the run:
1274	c "stats") only possible in 3D runs !
1275
1276	IF (callstats) THEN
1277
1278	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1279	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1280	call wstats(ngrid,"co2ice","CO2 ice cover",
1281	& "kg.m-2",2,co2ice)
1282	call wstats(ngrid,"fluxsurf_lw",
1283	& "Thermal IR radiative flux to surface","W.m-2",2,
1284	& fluxsurf_lw)
1285	call wstats(ngrid,"fluxsurf_sw",
1286	& "Solar radiative flux to surface","W.m-2",2,
1287	& fluxsurf_sw_tot)
1288	call wstats(ngrid,"fluxtop_lw",
1289	& "Thermal IR radiative flux to space","W.m-2",2,
1290	& fluxtop_lw)
1291	call wstats(ngrid,"fluxtop_sw",
1292	& "Solar radiative flux to space","W.m-2",2,
1293	& fluxtop_sw_tot)
1294	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1295	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1296	call wstats(ngrid,"v","Meridional (North-South) wind",
1297	& "m.s-1",3,zv)
1298	call wstats(ngrid,"w","Vertical (down-up) wind",
1299	& "m.s-1",3,pw)
1300	call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
1301	call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
1302	c call wstats(ngrid,"q2",
1303	c & "Boundary layer eddy kinetic energy",
1304	c & "m2.s-2",3,q2)
1305	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1306	c & emis)
1307	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
1308	c & 2,zstress)
1309	c call wstats(ngrid,"sw_htrt","sw heat.rate",
1310	c & "W.m-2",3,zdtsw)
1311	c call wstats(ngrid,"lw_htrt","lw heat.rate",
1312	c & "W.m-2",3,zdtlw)
1313
1314	if (tracer) then
1315	if (water) then
1316	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1317	& *mugaz/mmol(igcm_h2o_vap)
1318	call wstats(ngrid,"vmr_h2ovapor",
1319	& "H2O vapor volume mixing ratio","mol/mol",
1320	& 3,vmr)
1321	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1322	& *mugaz/mmol(igcm_h2o_ice)
1323	call wstats(ngrid,"vmr_h2oice",
1324	& "H2O ice volume mixing ratio","mol/mol",
1325	& 3,vmr)
1326	call wstats(ngrid,"h2o_ice_s",
1327	& "surface h2o_ice","kg/m2",
1328	& 2,qsurf(1,igcm_h2o_ice))
1329
1330	call wstats(ngrid,"mtot",
1331	& "total mass of water vapor","kg/m2",
1332	& 2,mtot)
1333	call wstats(ngrid,"icetot",
1334	& "total mass of water ice","kg/m2",
1335	& 2,icetot)
1336	call wstats(ngrid,"reffice",
1337	& "Mean reff","m",
1338	& 2,rave)
1339	c call wstats(ngrid,"rice",
1340	c & "Ice particle size","m",
1341	c & 3,rice)
1342	c If activice is true, tauTES is computed in aeropacity.F;
1343	if (.not.activice) then
1344	call wstats(ngrid,"tauTESap",
1345	& "tau abs 825 cm-1","",
1346	& 2,tauTES)
1347	endif
1348
1349	endif ! of if (water)
1350
1351	if (thermochem.or.photochem) then
1352	do iq=1,nq
1353	if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
1354	. (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
1355	. (noms(iq).eq."h2").or.
1356	. (noms(iq).eq."o3")) then
1357	do l=1,nlayer
1358	do ig=1,ngrid
1359	vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
1360	end do
1361	end do
1362	call wstats(ngrid,"vmr_"//trim(noms(iq)),
1363	. "Volume mixing ratio","mol/mol",3,vmr)
1364	endif
1365	enddo
1366	endif ! of if (thermochem.or.photochem)
1367
1368	endif ! of if (tracer)
1369
1370	IF(lastcall) THEN
1371	write (,) "Writing stats..."
1372	call mkstats(ierr)
1373	ENDIF
1374
1375	ENDIF !if callstats
1376
1377	c (Store EOF for Mars Climate database software)
1378	IF (calleofdump) THEN
1379	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1380	ENDIF
1381
1382
1383	#ifdef MESOSCALE
1384	!!!
1385	!!! OUTPUT FIELDS
1386	!!!
1387	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
1388	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
1389	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
1390	IF (igcm_h2o_ice .ne. 0) THEN
1391	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
1392	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1393	. * mugaz / mmol(igcm_h2o_ice)
1394	ENDIF
1395	c AUTOMATICALLY GENERATED FROM REGISTRY
1396	#include "fill_save.inc"
1397	#else
1398
1399	c ==========================================================
1400	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
1401	c any variable for diagnostic (output with period
1402	c "ecritphy", set in "run.def")
1403	c ==========================================================
1404	c WRITEDIAGFI can ALSO be called from any other subroutines
1405	c for any variables !!
1406	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1407	& emis)
1408	! call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
1409	! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
1410	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1411	& tsurf)
1412	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1413	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
1414	& co2ice)
1415	c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
1416	c & "K",2,zt(1,7))
1417	c call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1418	c & fluxsurf_lw)
1419	c call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1420	c & fluxsurf_sw_tot)
1421	c call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1422	c & fluxtop_lw)
1423	c call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1424	c & fluxtop_sw_tot)
1425	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1426	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1427	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1428	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1429	c call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1430	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1431	! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1432	c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
1433	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
1434	c & zstress)
1435	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
1436	c & 'w.m-2',3,zdtsw)
1437	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
1438	c & 'w.m-2',3,zdtlw)
1439	c CALL WRITEDIAGFI(ngridmx,'tauTESap',
1440	c & 'tau abs 825 cm-1',
1441	c & '',2,tauTES)
1442
1443	!!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL
1444	call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
1445	& "K",3,tsoil)
1446	call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
1447	& "K",3,inertiedat)
1448	!!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL
1449
1450	c ----------------------------------------------------------
1451	c Outputs of the CO2 cycle
1452	c ----------------------------------------------------------
1453
1454	if (tracer.and.(igcm_co2.ne.0)) then
1455	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
1456	! & "kg/kg",2,zq(1,1,igcm_co2))
1457	! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
1458	! & "kg/kg",3,zq(1,1,igcm_co2))
1459
1460	! Compute co2 column
1461	call zerophys(ngrid,co2col)
1462	do l=1,nlayermx
1463	do ig=1,ngrid
1464	co2col(ig)=co2col(ig)+
1465	& zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
1466	enddo
1467	enddo
1468	c call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
1469	c & co2col)
1470	endif ! of if (tracer.and.(igcm_co2.ne.0))
1471
1472	c ----------------------------------------------------------
1473	c Outputs of the water cycle
1474	c ----------------------------------------------------------
1475	if (tracer) then
1476	if (water) then
1477
1478	!!!! waterice = q01, voir readmeteo.F90
1479	call WRITEDIAGFI(ngridmx,'q01',noms(igcm_h2o_ice),
1480	& 'kg/kg',3,
1481	& zq(1:ngridmx,1:nlayermx,igcm_h2o_ice))
1482	!!!! watervapor = q02, voir readmeteo.F90
1483	call WRITEDIAGFI(ngridmx,'q02',noms(igcm_h2o_vap),
1484	& 'kg/kg',3,
1485	& zq(1:ngridmx,1:nlayermx,igcm_h2o_vap))
1486	!!!! surface waterice qsurf02 (voir readmeteo)
1487	call WRITEDIAGFI(ngridmx,'qsurf02','surface tracer',
1488	& 'kg.m-2',2,
1489	& qsurf(1:ngridmx,igcm_h2o_ice))
1490
1491	CALL WRITEDIAGFI(ngridmx,'mtot',
1492	& 'total mass of water vapor',
1493	& 'kg/m2',2,mtot)
1494	CALL WRITEDIAGFI(ngridmx,'icetot',
1495	& 'total mass of water ice',
1496	& 'kg/m2',2,icetot)
1497	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1498	c & *mugaz/mmol(igcm_h2o_ice)
1499	c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
1500	c & 'mol/mol',3,vmr)
1501	CALL WRITEDIAGFI(ngridmx,'reffice',
1502	& 'Mean reff',
1503	& 'm',2,rave)
1504	c call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
1505	c & 'm',3,rice)
1506	c If activice is true, tauTES is computed in aeropacity.F;
1507	if (.not.activice) then
1508	CALL WRITEDIAGFI(ngridmx,'tauTESap',
1509	& 'tau abs 825 cm-1',
1510	& '',2,tauTES)
1511	endif
1512	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
1513	& 'surface h2o_ice',
1514	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
1515	endif !(water)
1516
1517
1518	if (water.and..not.photochem) then
1519	iq=nq
1520	c write(str2(1:2),'(i2.2)') iq
1521	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
1522	c & 'kg.m-2',2,zdqscloud(1,iq))
1523	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
1524	c & 'kg/kg',3,zdqchim(1,1,iq))
1525	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
1526	c & 'kg/kg',3,zdqdif(1,1,iq))
1527	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
1528	c & 'kg/kg',3,zdqadj(1,1,iq))
1529	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
1530	c & 'kg/kg',3,zdqc(1,1,iq))
1531	endif !(water.and..not.photochem)
1532	endif
1533
1534	c ----------------------------------------------------------
1535	c Outputs of the dust cycle
1536	c ----------------------------------------------------------
1537
1538	c call WRITEDIAGFI(ngridmx,'tauref',
1539	c & 'Dust ref opt depth','NU',2,tauref)
1540
1541	if (tracer.and.(dustbin.ne.0)) then
1542	c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
1543	if (doubleq) then
1544	cc call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
1545	cc & 'kg.m-2',2,qsurf(1,1))
1546	cc call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
1547	cc & 'N.m-2',2,qsurf(1,2))
1548	cc call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
1549	cc & 'kg.m-2.s-1',2,zdqsdev(1,1))
1550	cc call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
1551	cc & 'kg.m-2.s-1',2,zdqssed(1,1))
1552	c call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
1553	c & 'm',3,rdust*ref_r0)
1554	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
1555	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
1556	call WRITEDIAGFI(ngridmx,'dustN','Dust number',
1557	& 'part/kg',3,pq(1,1,igcm_dust_number))
1558	else
1559	do iq=1,dustbin
1560	write(str2(1:2),'(i2.2)') iq
1561	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
1562	& 'kg/kg',3,zq(1,1,iq))
1563	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
1564	& 'kg.m-2',2,qsurf(1,iq))
1565	end do
1566	endif ! (doubleq)
1567	c if (submicron) then
1568	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
1569	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
1570	c endif ! (submicron)
1571	end if ! (tracer.and.(dustbin.ne.0))
1572
1573	c ----------------------------------------------------------
1574	c Outputs of thermals
1575	c ----------------------------------------------------------
1576	if (calltherm) then
1577
1578	! call WRITEDIAGFI(ngrid,'dtke',
1579	! & 'tendance tke thermiques','m2/s2',
1580	! & 3,dtke_th)
1581	! call WRITEDIAGFI(ngrid,'d_u_ajs',
1582	! & 'tendance u thermiques','m/s',
1583	! & 3,pdu_th*ptimestep)
1584	! call WRITEDIAGFI(ngrid,'d_v_ajs',
1585	! & 'tendance v thermiques','m/s',
1586	! & 3,pdv_th*ptimestep)
1587	! if (tracer) then
1588	! if (nq .eq. 2) then
1589	! call WRITEDIAGFI(ngrid,'deltaq_th',
1590	! & 'delta q thermiques','kg/kg',
1591	! & 3,ptimestep*pdq_th(:,:,2))
1592	! endif
1593	! endif
1594
1595	lmax_th_out(:)=real(lmax_th(:))
1596
1597	call WRITEDIAGFI(ngridmx,'lmax_th',
1598	& 'hauteur du thermique','K',
1599	& 2,lmax_th_out)
1600	call WRITEDIAGFI(ngridmx,'hfmax_th',
1601	& 'maximum TH heat flux','K.m/s',
1602	& 2,hfmax_th)
1603	call WRITEDIAGFI(ngridmx,'wmax_th',
1604	& 'maximum TH vertical velocity','m/s',
1605	& 2,wmax_th)
1606
1607	endif
1608
1609	c ----------------------------------------------------------
1610	c Output in netcdf file "diagsoil.nc" for subterranean
1611	c variables (output every "ecritphy", as for writediagfi)
1612	c ----------------------------------------------------------
1613
1614	! Write soil temperature
1615	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
1616	! & 3,tsoil)
1617	! Write surface temperature
1618	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
1619	! & 2,tsurf)
1620
1621	c ==========================================================
1622	c END OF WRITEDIAGFI
1623	c ==========================================================
1624	#endif
1625
1626	ELSE ! if(ngrid.eq.1)
1627
1628	print,'Ls =',zls180./pi,
1629	& ' tauref(700 Pa) =',tauref
1630	c ----------------------------------------------------------------------
1631	c Output in grads file "g1d" (ONLY when using testphys1d)
1632	c (output at every X physical timestep)
1633	c ----------------------------------------------------------------------
1634	c
1635	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
1636	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
1637	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
1638
1639	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
1640	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
1641	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
1642	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
1643
1644	! THERMALS STUFF 1D
1645
1646	if(calltherm) then
1647
1648	lmax_th_out(:)=real(lmax_th(:))
1649
1650	call WRITEDIAGFI(ngridmx,'lmax_th',
1651	& 'hauteur du thermique','point',
1652	& 0,lmax_th_out)
1653	call WRITEDIAGFI(ngridmx,'hfmax_th',
1654	& 'maximum TH heat flux','K.m/s',
1655	& 0,hfmax_th)
1656	call WRITEDIAGFI(ngridmx,'wmax_th',
1657	& 'maximum TH vertical velocity','m/s',
1658	& 0,wmax_th)
1659
1660
1661	co2col(:)=0.
1662	dummycol(:)=0.
1663	if (tracer) then
1664	do l=1,nlayermx
1665	do ig=1,ngrid
1666	co2col(ig)=co2col(ig)+
1667	& zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
1668	if (nqmx .gt. 1) then
1669	iq=2 ! to avoid out-of-bounds spotting by picky compilers
1670	! when gcm is compiled with only one tracer
1671	dummycol(ig)=dummycol(ig)+
1672	& zq(ig,l,iq)*(pplev(ig,l)-pplev(ig,l+1))/g
1673	endif
1674	enddo
1675	enddo
1676
1677	end if
1678	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
1679	& ,'kg/m-2',0,co2col)
1680	call WRITEDIAGFI(ngrid,'dummycol','integrated dummy mass' &
1681	& ,'kg/m-2',0,dummycol)
1682	endif
1683	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
1684	& ,'m/s',1,pw)
1685	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
1686	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
1687	& tsurf)
1688
1689	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
1690	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
1691	! or output in diagfi.nc (for testphys1d)
1692	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
1693	call WRITEDIAGFI(ngridmx,'temp','Temperature',
1694	& 'K',1,zt)
1695
1696	if(tracer) then
1697	c CALL writeg1d(ngrid,1,tau,'tau','SI')
1698	do iq=1,nq
1699	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
1700	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
1701	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
1702	end do
1703	end if
1704
1705	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
1706
1707	do l=2,nlayer-1
1708	tmean=zt(1,l)
1709	if(zt(1,l).ne.zt(1,l-1))
1710	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
1711	zlocal(l)= zlocal(l-1)
1712	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
1713	enddo
1714	zlocal(nlayer)= zlocal(nlayer-1)-
1715	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
1716	& rnew(1,nlayer)*tmean/g
1717
1718	END IF ! if(ngrid.ne.1)
1719
1720	icount=icount+1
1721	RETURN
1722	END

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

Download in other formats:

Original Format