Context Navigation

physiq.F @ 283

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

LMDZ.MARS:

AS: J'ai fait ce commit et fait a la main un merge avec les modifications entre temps

24/08/11 == TN

Attempts to tune the water cycle by adding outliers

+ A few structural changes !!

watercap.h is now obsolete and removed -- all is in surfdat.h
watercaptag initialized in surfini.F (up to 7 areas defined) instead of initracer.F
- settings proposed by AS commented
- experiments by TN decommented. use with caution.
water ice albedo and thermal inertia in callphys.def and inifis.F
water ice albedo in surfini.F
water ice albedo computation in albedocaps.F90
alb_surfice is now obsolete in physiq.F, albedo_h2o_ice is used instead
frost_albedo_threshold defined in surfdat.h
water ice thermal inertia in soil.F

TODO: * calibrate thermal inertia and ice albedo

have a look at subgrid-scale ice with dryness ?

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

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

Download in other formats:

Original Format