Context Navigation

physiq.F @ 324

Last change on this file since 324 was 319, checked in by aslmd, 14 years ago
LMDZ.MARS : minor changes to thermals [output L_mo + default settings]. transparent to users.
File size: 66.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 "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 [meters]
330	REAL ustar(ngridmx),tstar(ngridmx) ! friction velocity and friction potential temp
331	REAL L_mo(ngridmx)
332	REAL zu2(ngridmx)
333	c=======================================================================
334
335	c 1. Initialisation:
336	c -----------------
337
338	c 1.1 Initialisation only at first call
339	c ---------------------------------------
340	IF (firstcall) THEN
341
342	c variables set to 0
343	c ~~~~~~~~~~~~~~~~~~
344	aerosol(:,:,:)=0
345	dtrad(:,:)=0
346	fluxrad(:)=0
347
348	wmax_th(:)=0.
349
350	c read startfi
351	c ~~~~~~~~~~~~
352	#ifndef MESOSCALE
353	! Read netcdf initial physical parameters.
354	CALL phyetat0 ("startfi.nc",0,0,
355	& nsoilmx,nq,
356	& day_ini,time_phys,
357	& tsurf,tsoil,emis,q2,qsurf,co2ice)
358	#else
359	#include "meso_inc/meso_inc_ini.F"
360	#endif
361
362	if (pday.ne.day_ini) then
363	write(,) "PHYSIQ: ERROR: bad synchronization between ",
364	& "physics and dynamics"
365	write(,) "dynamics day: ",pday
366	write(,) "physics day: ",day_ini
367	stop
368	endif
369
370	write (,) 'In physiq day_ini =', day_ini
371
372	c initialize tracers
373	c ~~~~~~~~~~~~~~~~~~
374	tracerdyn=tracer
375	IF (tracer) THEN
376	CALL initracer(qsurf,co2ice)
377	ENDIF ! end tracer
378
379	c Initialize albedo and orbital calculation
380	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
381	CALL surfini(ngrid,co2ice,qsurf,albedo)
382	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
383
384	c initialize soil
385	c ~~~~~~~~~~~~~~~
386	IF (callsoil) THEN
387	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
388	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
389	ELSE
390	PRINT*,
391	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
392	DO ig=1,ngrid
393	capcal(ig)=1.e5
394	fluxgrd(ig)=0.
395	ENDDO
396	ENDIF
397	icount=1
398
399	#ifndef MESOSCALE
400	c Initialize thermospheric parameters
401	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
402
403	if (callthermos) call param_read
404	#endif
405	c Initialize R and Cp as constant
406
407	if (.not.callthermos .and. .not.photochem) then
408	do l=1,nlayermx
409	do ig=1,ngridmx
410	rnew(ig,l)=r
411	cpnew(ig,l)=cpp
412	mmean(ig,l)=mugaz
413	enddo
414	enddo
415	endif
416
417	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
418	write(,)"physiq: water_param Surface water ice albedo:",
419	. albedo_h2o_ice
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	pdv(:,:)=0
439	pdu(:,:)=0
440	pdt(:,:)=0
441	pdq(:,:,:)=0
442	pdpsrf(:)=0
443	zflubid(:)=0
444	zdtsurf(:)=0
445	dqsurf(:,:)=0
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, 'Ls =',zls180./pi,
552	& 'check lat lon', lati(igout)*180/pi,
553	& long(igout)*180/pi
554	print*, 'tauref(700 Pa) =',tauref(igout),
555	& ' 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	zdtnirco2(:,:)=0
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	zdum1(:,:)=0
685	zdum2(:,:)=0
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
693	#ifdef MESOSCALE
694	IF (.not.flag_LES) THEN
695	#endif
696	c ----------------------
697	c Treatment of a special case : using new surface layer (Richardson based)
698	c without using the thermals in gcm and mesoscale can yield problems in
699	c weakly unstable situations when winds are near to 0. For those cases, we add
700	c a unit subgrid gustiness. Remember that thermals should be used we using the
701	c Richardson based surface layer model.
702	IF ( .not.calltherm .and. callrichsl ) THEN
703	DO ig=1, ngridmx
704	IF (zh(ig,1) .lt. tsurf(ig)) THEN
705	wmax_th(ig)=1.
706	ENDIF
707	ENDDO
708	ENDIF
709	c ----------------------
710	#ifdef MESOSCALE
711	ENDIF
712	#endif
713
714
715	c Calling vdif (Martian version WITH CO2 condensation)
716	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
717	$ ptimestep,capcal,lwrite,
718	$ pplay,pplev,zzlay,zzlev,z0,
719	$ pu,pv,zh,pq,tsurf,emis,qsurf,
720	$ zdum1,zdum2,zdh,pdq,zflubid,
721	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
722	& zdqdif,zdqsdif,wmax_th,zcdv,zcdh)
723
724	#ifdef MESOSCALE
725	#include "meso_inc/meso_inc_les.F"
726	#endif
727	DO l=1,nlayer
728	DO ig=1,ngrid
729	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
730	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
731	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
732
733	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
734
735	ENDDO
736	ENDDO
737
738	DO ig=1,ngrid
739	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
740	ENDDO
741
742	if (tracer) then
743	DO iq=1, nq
744	DO l=1,nlayer
745	DO ig=1,ngrid
746	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
747	ENDDO
748	ENDDO
749	ENDDO
750	DO iq=1, nq
751	DO ig=1,ngrid
752	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
753	ENDDO
754	ENDDO
755	end if ! of if (tracer)
756
757	ELSE
758	DO ig=1,ngrid
759	zdtsurf(ig)=zdtsurf(ig)+
760	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
761	ENDDO
762	#ifdef MESOSCALE
763	IF (flag_LES) THEN
764	write(,) 'LES mode !'
765	write(,) 'Please set calldifv to T in callphys.def'
766	STOP
767	ENDIF
768	#endif
769	ENDIF ! of IF (calldifv)
770
771	c-----------------------------------------------------------------------
772	c TEST. Thermals :
773	c HIGHLY EXPERIMENTAL, BEWARE !!
774	c -----------------------------
775
776	if(calltherm) then
777
778	call calltherm_interface(firstcall,
779	$ long,lati,zzlev,zzlay,
780	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
781	$ pplay,pplev,pphi,zpopsk,
782	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
783	$ dtke_th,hfmax_th,wmax_th)
784
785	DO l=1,nlayer
786	DO ig=1,ngrid
787	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
788	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
789	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
790	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
791	ENDDO
792	ENDDO
793
794	DO ig=1,ngrid
795	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
796	ENDDO
797
798	if (tracer) then
799	DO iq=1,nq
800	DO l=1,nlayer
801	DO ig=1,ngrid
802	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
803	ENDDO
804	ENDDO
805	ENDDO
806	endif
807
808	lmax_th_out(:)=real(lmax_th(:))
809
810	else !of if calltherm
811	lmax_th(:)=0
812	wmax_th(:)=0.
813	lmax_th_out(:)=0.
814	end if
815
816	c-----------------------------------------------------------------------
817	c 5. Dry convective adjustment:
818	c -----------------------------
819
820	IF(calladj) THEN
821
822	DO l=1,nlayer
823	DO ig=1,ngrid
824	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
825	ENDDO
826	ENDDO
827	zduadj(:,:)=0
828	zdvadj(:,:)=0
829	zdhadj(:,:)=0
830
831	CALL convadj(ngrid,nlayer,nq,ptimestep,
832	$ pplay,pplev,zpopsk,lmax_th,
833	$ pu,pv,zh,pq,
834	$ pdu,pdv,zdh,pdq,
835	$ zduadj,zdvadj,zdhadj,
836	$ zdqadj)
837
838
839	DO l=1,nlayer
840	DO ig=1,ngrid
841	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
842	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
843	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
844
845	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
846	ENDDO
847	ENDDO
848
849	if(tracer) then
850	DO iq=1, nq
851	DO l=1,nlayer
852	DO ig=1,ngrid
853	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
854	ENDDO
855	ENDDO
856	ENDDO
857	end if
858	ENDIF ! of IF(calladj)
859
860	c-----------------------------------------------------------------------
861	c 6. Carbon dioxide condensation-sublimation:
862	c -------------------------------------------
863
864	IF (callcond) THEN
865	CALL newcondens(ngrid,nlayer,nq,ptimestep,
866	$ capcal,pplay,pplev,tsurf,pt,
867	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
868	$ co2ice,albedo,emis,
869	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
870	$ fluxsurf_sw,zls)
871
872	DO l=1,nlayer
873	DO ig=1,ngrid
874	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
875	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
876	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
877	ENDDO
878	ENDDO
879	DO ig=1,ngrid
880	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
881	ENDDO
882
883	IF (tracer) THEN
884	DO iq=1, nq
885	DO l=1,nlayer
886	DO ig=1,ngrid
887	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
888	ENDDO
889	ENDDO
890	ENDDO
891	ENDIF ! of IF (tracer)
892
893	ENDIF ! of IF (callcond)
894
895	c-----------------------------------------------------------------------
896	c 7. Specific parameterizations for tracers
897	c: -----------------------------------------
898
899	if (tracer) then
900
901	c 7a. Water and ice
902	c ---------------
903
904	c ---------------------------------------
905	c Water ice condensation in the atmosphere
906	c ----------------------------------------
907	IF (water) THEN
908
909	call watercloud(ngrid,nlayer,ptimestep,
910	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
911	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
912	& nq,naerkind,tau,
913	& ccn,rdust,rice,nuice)
914	if (activice) then
915	c Temperature variation due to latent heat release
916	DO l=1,nlayer
917	DO ig=1,ngrid
918	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
919	ENDDO
920	ENDDO
921	endif
922
923	! increment water vapour and ice atmospheric tracers tendencies
924	IF (water) THEN
925	DO l=1,nlayer
926	DO ig=1,ngrid
927	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
928	& zdqcloud(ig,l,igcm_h2o_vap)
929	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
930	& zdqcloud(ig,l,igcm_h2o_ice)
931	ENDDO
932	ENDDO
933	ENDIF ! of IF (water) THEN
934	! Increment water ice surface tracer tendency
935	DO ig=1,ngrid
936	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
937	& zdqscloud(ig,igcm_h2o_ice)
938	ENDDO
939
940	END IF ! of IF (water)
941
942	c 7b. Chemical species
943	c ------------------
944
945	#ifndef MESOSCALE
946	c --------------
947	c photochemistry :
948	c --------------
949	IF (photochem .or. thermochem) then
950	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
951	& zzlay,zday,pq,pdq,rice,
952	& zdqchim,zdqschim,zdqcloud,zdqscloud)
953	!NB: Photochemistry includes condensation of H2O2
954
955	! increment values of tracers:
956	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
957	! tracers is zero anyways
958	DO l=1,nlayer
959	DO ig=1,ngrid
960	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
961	ENDDO
962	ENDDO
963	ENDDO ! of DO iq=1,nq
964	! add condensation tendency for H2O2
965	if (igcm_h2o2.ne.0) then
966	DO l=1,nlayer
967	DO ig=1,ngrid
968	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
969	& +zdqcloud(ig,l,igcm_h2o2)
970	ENDDO
971	ENDDO
972	endif
973
974	! increment surface values of tracers:
975	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
976	! tracers is zero anyways
977	DO ig=1,ngrid
978	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
979	ENDDO
980	ENDDO ! of DO iq=1,nq
981	! add condensation tendency for H2O2
982	if (igcm_h2o2.ne.0) then
983	DO ig=1,ngrid
984	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
985	& +zdqscloud(ig,igcm_h2o2)
986	ENDDO
987	endif
988
989	END IF ! of IF (photochem.or.thermochem)
990	#endif
991
992	c 7c. Aerosol particles
993	c -------------------
994
995	c ----------
996	c Dust devil :
997	c ----------
998	IF(callddevil) then
999	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
1000	& zdqdev,zdqsdev)
1001
1002	if (dustbin.ge.1) then
1003	do iq=1,nq
1004	DO l=1,nlayer
1005	DO ig=1,ngrid
1006	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
1007	ENDDO
1008	ENDDO
1009	enddo
1010	do iq=1,nq
1011	DO ig=1,ngrid
1012	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
1013	ENDDO
1014	enddo
1015	endif ! of if (dustbin.ge.1)
1016
1017	END IF ! of IF (callddevil)
1018
1019	c -------------
1020	c Sedimentation : acts also on water ice
1021	c -------------
1022	IF (sedimentation) THEN
1023	!call zerophys(ngridnlayernq, zdqsed)
1024	zdqsed(1:ngrid,1:nlayer,1:nq)=0
1025	!call zerophys(ngrid*nq, zdqssed)
1026	zdqssed(1:ngrid,1:nq)=0
1027
1028	call callsedim(ngrid,nlayer, ptimestep,
1029	& pplev,zzlev, pt, rdust, rice,
1030	& pq, pdq, zdqsed, zdqssed,nq)
1031
1032	DO iq=1, nq
1033	DO l=1,nlayer
1034	DO ig=1,ngrid
1035	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1036	ENDDO
1037	ENDDO
1038	ENDDO
1039	DO iq=1, nq
1040	DO ig=1,ngrid
1041	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1042	ENDDO
1043	ENDDO
1044	END IF ! of IF (sedimentation)
1045
1046	c 7d. Updates
1047	c ---------
1048
1049	DO iq=1, nq
1050	DO ig=1,ngrid
1051
1052	c ---------------------------------
1053	c Updating tracer budget on surface
1054	c ---------------------------------
1055	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1056
1057	ENDDO ! (ig)
1058	ENDDO ! (iq)
1059
1060	endif ! of if (tracer)
1061
1062	#ifndef MESOSCALE
1063	c-----------------------------------------------------------------------
1064	c 8. THERMOSPHERE CALCULATION
1065	c-----------------------------------------------------------------------
1066
1067	if (callthermos) then
1068	call thermosphere(pplev,pplay,dist_sol,
1069	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1070	& pt,pq,pu,pv,pdt,pdq,
1071	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1072
1073	DO l=1,nlayer
1074	DO ig=1,ngrid
1075	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1076	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1077	& +zdteuv(ig,l)
1078	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1079	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1080	DO iq=1, nq
1081	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1082	ENDDO
1083	ENDDO
1084	ENDDO
1085
1086	endif ! of if (callthermos)
1087	#endif
1088
1089	c-----------------------------------------------------------------------
1090	c 9. Surface and sub-surface soil temperature
1091	c-----------------------------------------------------------------------
1092	c
1093	c
1094	c 9.1 Increment Surface temperature:
1095	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1096
1097	DO ig=1,ngrid
1098	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1099	ENDDO
1100
1101	c Prescribe a cold trap at south pole (except at high obliquity !!)
1102	c Temperature at the surface is set there to be the temperature
1103	c corresponding to equilibrium temperature between phases of CO2
1104
1105	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
1106	#ifndef MESOSCALE
1107	if (caps.and.(obliquit.lt.27.)) then
1108	! NB: Updated surface pressure, at grid point 'ngrid', is
1109	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1110	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1111	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1112	endif
1113	#endif
1114	c -------------------------------------------------------------
1115	c Change of surface albedo in case of ground frost
1116	c everywhere except on the north permanent cap and in regions
1117	c covered by dry ice.
1118	c ALWAYS PLACE these lines after newcondens !!!
1119	c -------------------------------------------------------------
1120	do ig=1,ngrid
1121	if ((co2ice(ig).eq.0).and.
1122	& (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
1123	albedo(ig,1) = albedo_h2o_ice
1124	albedo(ig,2) = albedo_h2o_ice
1125	c write(,) "frost thickness", qsurf(ig,igcm_h2o_ice)
1126	c write(,) "physiq.F frost :"
1127	c & ,lati(ig)180./pi, long(ig)180./pi
1128	endif
1129	enddo ! of do ig=1,ngrid
1130	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
1131
1132	c
1133	c 9.2 Compute soil temperatures and subsurface heat flux:
1134	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1135	IF (callsoil) THEN
1136	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1137	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
1138	ENDIF
1139
1140	c-----------------------------------------------------------------------
1141	c 10. Write output files
1142	c ----------------------
1143
1144	c -------------------------------
1145	c Dynamical fields incrementation
1146	c -------------------------------
1147	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1148	! temperature, zonal and meridional wind
1149	DO l=1,nlayer
1150	DO ig=1,ngrid
1151	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1152	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1153	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1154	ENDDO
1155	ENDDO
1156
1157	! tracers
1158	DO iq=1, nq
1159	DO l=1,nlayer
1160	DO ig=1,ngrid
1161	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1162	ENDDO
1163	ENDDO
1164	ENDDO
1165
1166	! surface pressure
1167	DO ig=1,ngrid
1168	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1169	ENDDO
1170
1171	! pressure
1172	DO l=1,nlayer
1173	DO ig=1,ngrid
1174	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1175	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1176	ENDDO
1177	ENDDO
1178
1179	! Density
1180	DO l=1,nlayer
1181	DO ig=1,ngrid
1182	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1183	ENDDO
1184	ENDDO
1185
1186	! Potential Temperature
1187
1188	DO ig=1,ngridmx
1189	DO l=1,nlayermx
1190	zh(ig,l) = zt(ig,l)(zplay(ig,l)/zplev(ig,1))*rcp
1191	ENDDO
1192	ENDDO
1193
1194
1195	c Compute surface stress : (NB: z0 is a common in surfdat.h)
1196	c DO ig=1,ngrid
1197	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
1198	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1199	c ENDDO
1200
1201	c Sum of fluxes in solar spectral bands (for output only)
1202	DO ig=1,ngrid
1203	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1204	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1205	ENDDO
1206	c ***** TEST ****************************************************
1207	ztim1 = 999
1208	DO l=1,nlayer
1209	DO ig=1,ngrid
1210	if (pt(ig,l).lt.ztim1) then
1211	ztim1 = pt(ig,l)
1212	igmin = ig
1213	lmin = l
1214	end if
1215	ENDDO
1216	ENDDO
1217	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1218	write(,) 'PHYSIQ: stability WARNING :'
1219	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1220	& 'ig l =', igmin, lmin
1221	end if
1222	c *******************************************************************
1223
1224	c ---------------------
1225	c Outputs to the screen
1226	c ---------------------
1227
1228	IF (lwrite) THEN
1229	PRINT*,'Global diagnostics for the physics'
1230	PRINT,'Variables and their increments x and dx/dt dt'
1231	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1232	WRITE(*,'(2f10.5,2f15.5)')
1233	s tsurf(igout),zdtsurf(igout)*ptimestep,
1234	s pplev(igout,1),pdpsrf(igout)*ptimestep
1235	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1236	WRITE(*,'(i4,6f10.5)') (l,
1237	s pu(igout,l),pdu(igout,l)*ptimestep,
1238	s pv(igout,l),pdv(igout,l)*ptimestep,
1239	s pt(igout,l),pdt(igout,l)*ptimestep,
1240	s l=1,nlayer)
1241	ENDIF ! of IF (lwrite)
1242
1243	IF (ngrid.NE.1) THEN
1244
1245	#ifndef MESOSCALE
1246	c -------------------------------------------------------------------
1247	c Writing NetCDF file "RESTARTFI" at the end of the run
1248	c -------------------------------------------------------------------
1249	c Note: 'restartfi' is stored just before dynamics are stored
1250	c in 'restart'. Between now and the writting of 'restart',
1251	c there will have been the itau=itau+1 instruction and
1252	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1253	c thus we store for time=time+dtvr
1254
1255	IF(lastcall) THEN
1256	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1257	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1258	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1259	. ptimestep,pday,
1260	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
1261	. area,albedodat,inertiedat,zmea,zstd,zsig,
1262	. zgam,zthe)
1263	ENDIF
1264	#endif
1265
1266	c -------------------------------------------------------------------
1267	c Calculation of diagnostic variables written in both stats and
1268	c diagfi files
1269	c -------------------------------------------------------------------
1270
1271	if (tracer) then
1272	if (water) then
1273
1274	mtot(:)=0
1275	icetot(:)=0
1276	rave(:)=0
1277	tauTES(:)=0
1278	do ig=1,ngrid
1279	do l=1,nlayermx
1280	mtot(ig) = mtot(ig) +
1281	& zq(ig,l,igcm_h2o_vap) *
1282	& (pplev(ig,l) - pplev(ig,l+1)) / g
1283	icetot(ig) = icetot(ig) +
1284	& zq(ig,l,igcm_h2o_ice) *
1285	& (pplev(ig,l) - pplev(ig,l+1)) / g
1286	rave(ig) = rave(ig) +
1287	& zq(ig,l,igcm_h2o_ice) *
1288	& (pplev(ig,l) - pplev(ig,l+1)) / g *
1289	& rice(ig,l) * (1.+nuice_ref)
1290	c Computing abs optical depth at 825 cm-1 in each
1291	c layer to simulate NEW TES retrieval
1292	Qabsice = min(
1293	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
1294	& )
1295	opTES(ig,l)= 0.75 * Qabsice *
1296	& zq(ig,l,igcm_h2o_ice) *
1297	& (pplev(ig,l) - pplev(ig,l+1)) / g
1298	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1299	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1300	enddo
1301	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
1302	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1303	enddo
1304
1305	endif ! of if (water)
1306	endif ! of if (tracer)
1307
1308	c -----------------------------------------------------------------
1309	c WSTATS: Saving statistics
1310	c -----------------------------------------------------------------
1311	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1312	c which can later be used to make the statistic files of the run:
1313	c "stats") only possible in 3D runs !
1314
1315	IF (callstats) THEN
1316
1317	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1318	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1319	call wstats(ngrid,"co2ice","CO2 ice cover",
1320	& "kg.m-2",2,co2ice)
1321	call wstats(ngrid,"fluxsurf_lw",
1322	& "Thermal IR radiative flux to surface","W.m-2",2,
1323	& fluxsurf_lw)
1324	call wstats(ngrid,"fluxsurf_sw",
1325	& "Solar radiative flux to surface","W.m-2",2,
1326	& fluxsurf_sw_tot)
1327	call wstats(ngrid,"fluxtop_lw",
1328	& "Thermal IR radiative flux to space","W.m-2",2,
1329	& fluxtop_lw)
1330	call wstats(ngrid,"fluxtop_sw",
1331	& "Solar radiative flux to space","W.m-2",2,
1332	& fluxtop_sw_tot)
1333	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1334	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1335	call wstats(ngrid,"v","Meridional (North-South) wind",
1336	& "m.s-1",3,zv)
1337	call wstats(ngrid,"w","Vertical (down-up) wind",
1338	& "m.s-1",3,pw)
1339	call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
1340	call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
1341	c call wstats(ngrid,"q2",
1342	c & "Boundary layer eddy kinetic energy",
1343	c & "m2.s-2",3,q2)
1344	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1345	c & emis)
1346	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
1347	c & 2,zstress)
1348	c call wstats(ngrid,"sw_htrt","sw heat.rate",
1349	c & "W.m-2",3,zdtsw)
1350	c call wstats(ngrid,"lw_htrt","lw heat.rate",
1351	c & "W.m-2",3,zdtlw)
1352
1353	if (tracer) then
1354	if (water) then
1355	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1356	& *mugaz/mmol(igcm_h2o_vap)
1357	call wstats(ngrid,"vmr_h2ovapor",
1358	& "H2O vapor volume mixing ratio","mol/mol",
1359	& 3,vmr)
1360	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1361	& *mugaz/mmol(igcm_h2o_ice)
1362	call wstats(ngrid,"vmr_h2oice",
1363	& "H2O ice volume mixing ratio","mol/mol",
1364	& 3,vmr)
1365	call wstats(ngrid,"h2o_ice_s",
1366	& "surface h2o_ice","kg/m2",
1367	& 2,qsurf(1,igcm_h2o_ice))
1368
1369	call wstats(ngrid,"mtot",
1370	& "total mass of water vapor","kg/m2",
1371	& 2,mtot)
1372	call wstats(ngrid,"icetot",
1373	& "total mass of water ice","kg/m2",
1374	& 2,icetot)
1375	call wstats(ngrid,"reffice",
1376	& "Mean reff","m",
1377	& 2,rave)
1378	c call wstats(ngrid,"rice",
1379	c & "Ice particle size","m",
1380	c & 3,rice)
1381	c If activice is true, tauTES is computed in aeropacity.F;
1382	if (.not.activice) then
1383	call wstats(ngrid,"tauTESap",
1384	& "tau abs 825 cm-1","",
1385	& 2,tauTES)
1386	endif
1387
1388	endif ! of if (water)
1389
1390	if (thermochem.or.photochem) then
1391	do iq=1,nq
1392	if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
1393	. (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
1394	. (noms(iq).eq."h2").or.
1395	. (noms(iq).eq."o3")) then
1396	do l=1,nlayer
1397	do ig=1,ngrid
1398	vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
1399	end do
1400	end do
1401	call wstats(ngrid,"vmr_"//trim(noms(iq)),
1402	. "Volume mixing ratio","mol/mol",3,vmr)
1403	endif
1404	enddo
1405	endif ! of if (thermochem.or.photochem)
1406
1407	endif ! of if (tracer)
1408
1409	IF(lastcall) THEN
1410	write (,) "Writing stats..."
1411	call mkstats(ierr)
1412	ENDIF
1413
1414	ENDIF !if callstats
1415
1416	c (Store EOF for Mars Climate database software)
1417	IF (calleofdump) THEN
1418	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1419	ENDIF
1420
1421
1422	#ifdef MESOSCALE
1423	!!!
1424	!!! OUTPUT FIELDS
1425	!!!
1426	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
1427	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
1428	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
1429	!! JF
1430	TAU_lay(:)=tau(:,1)!!true opacity (not a reference like tauref)
1431	IF (igcm_dust_mass .ne. 0) THEN
1432	qsurfice_dust(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
1433	ENDIF
1434	IF (igcm_h2o_ice .ne. 0) THEN
1435	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
1436	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1437	. * mugaz / mmol(igcm_h2o_ice)
1438	ENDIF
1439	!! Dust quantity integration along the vertical axe
1440	dustot(:)=0
1441	do ig=1,ngrid
1442	do l=1,nlayermx
1443	dustot(ig) = dustot(ig) +
1444	& zq(ig,l,igcm_dust_mass)
1445	& * (pplev(ig,l) - pplev(ig,l+1)) / g
1446	enddo
1447	enddo
1448	c AUTOMATICALLY GENERATED FROM REGISTRY
1449	#include "fill_save.inc"
1450	#else
1451
1452	c ==========================================================
1453	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
1454	c any variable for diagnostic (output with period
1455	c "ecritphy", set in "run.def")
1456	c ==========================================================
1457	c WRITEDIAGFI can ALSO be called from any other subroutines
1458	c for any variables !!
1459	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1460	c & emis)
1461	! call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
1462	! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
1463	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1464	& tsurf)
1465	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1466	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
1467	& co2ice)
1468
1469	c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
1470	c & "K",2,zt(1,7))
1471	c call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1472	c & fluxsurf_lw)
1473	c call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1474	c & fluxsurf_sw_tot)
1475	c call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1476	c & fluxtop_lw)
1477	c call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1478	c & fluxtop_sw_tot)
1479	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1480	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1481	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1482	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1483	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1484	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1485	! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1486	c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
1487	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
1488	c & zstress)
1489	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
1490	c & 'w.m-2',3,zdtsw)
1491	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
1492	c & 'w.m-2',3,zdtlw)
1493	c CALL WRITEDIAGFI(ngridmx,'tauTESap',
1494	c & 'tau abs 825 cm-1',
1495	c & '',2,tauTES)
1496
1497	#ifdef MESOINI
1498	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1499	& emis)
1500	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1501	call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
1502	& "K",3,tsoil)
1503	call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
1504	& "K",3,inertiedat)
1505	#endif
1506
1507
1508	c ----------------------------------------------------------
1509	c Outputs of the CO2 cycle
1510	c ----------------------------------------------------------
1511
1512	if (tracer.and.(igcm_co2.ne.0)) then
1513	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
1514	! & "kg/kg",2,zq(1,1,igcm_co2))
1515	! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
1516	! & "kg/kg",3,zq(1,1,igcm_co2))
1517
1518	! Compute co2 column
1519	co2col(:)=0
1520	do l=1,nlayermx
1521	do ig=1,ngrid
1522	co2col(ig)=co2col(ig)+
1523	& zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
1524	enddo
1525	enddo
1526	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
1527	& co2col)
1528	endif ! of if (tracer.and.(igcm_co2.ne.0))
1529
1530	c ----------------------------------------------------------
1531	c Outputs of the water cycle
1532	c ----------------------------------------------------------
1533	if (tracer) then
1534	if (water) then
1535
1536	#ifdef MESOINI
1537	!!!! waterice = q01, voir readmeteo.F90
1538	call WRITEDIAGFI(ngridmx,'q01',noms(igcm_h2o_ice),
1539	& 'kg/kg',3,
1540	& zq(1:ngridmx,1:nlayermx,igcm_h2o_ice))
1541	!!!! watervapor = q02, voir readmeteo.F90
1542	call WRITEDIAGFI(ngridmx,'q02',noms(igcm_h2o_vap),
1543	& 'kg/kg',3,
1544	& zq(1:ngridmx,1:nlayermx,igcm_h2o_vap))
1545	!!!! surface waterice qsurf02 (voir readmeteo)
1546	call WRITEDIAGFI(ngridmx,'qsurf02','surface tracer',
1547	& 'kg.m-2',2,
1548	& qsurf(1:ngridmx,igcm_h2o_ice))
1549	#endif
1550
1551	CALL WRITEDIAGFI(ngridmx,'mtot',
1552	& 'total mass of water vapor',
1553	& 'kg/m2',2,mtot)
1554	CALL WRITEDIAGFI(ngridmx,'icetot',
1555	& 'total mass of water ice',
1556	& 'kg/m2',2,icetot)
1557	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1558	c & *mugaz/mmol(igcm_h2o_ice)
1559	c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
1560	c & 'mol/mol',3,vmr)
1561	CALL WRITEDIAGFI(ngridmx,'reffice',
1562	& 'Mean reff',
1563	& 'm',2,rave)
1564	c call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
1565	c & 'm',3,rice)
1566	c If activice is true, tauTES is computed in aeropacity.F;
1567	if (.not.activice) then
1568	CALL WRITEDIAGFI(ngridmx,'tauTESap',
1569	& 'tau abs 825 cm-1',
1570	& '',2,tauTES)
1571	endif
1572	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
1573	& 'surface h2o_ice',
1574	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
1575
1576	if (caps) then
1577	do ig=1,ngridmx
1578	if (watercaptag(ig)) watercapflag(ig) = 1
1579	enddo
1580	CALL WRITEDIAGFI(ngridmx,'watercaptag',
1581	& 'Ice water caps',
1582	& '',2,watercapflag)
1583	endif
1584	CALL WRITEDIAGFI(ngridmx,'albedo',
1585	& 'albedo',
1586	& '',2,albedo(1:ngridmx,1))
1587	endif !(water)
1588
1589
1590	if (water.and..not.photochem) then
1591	iq=nq
1592	c write(str2(1:2),'(i2.2)') iq
1593	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
1594	c & 'kg.m-2',2,zdqscloud(1,iq))
1595	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
1596	c & 'kg/kg',3,zdqchim(1,1,iq))
1597	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
1598	c & 'kg/kg',3,zdqdif(1,1,iq))
1599	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
1600	c & 'kg/kg',3,zdqadj(1,1,iq))
1601	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
1602	c & 'kg/kg',3,zdqc(1,1,iq))
1603	endif !(water.and..not.photochem)
1604	endif
1605
1606	c ----------------------------------------------------------
1607	c Outputs of the dust cycle
1608	c ----------------------------------------------------------
1609
1610	c call WRITEDIAGFI(ngridmx,'tauref',
1611	c & 'Dust ref opt depth','NU',2,tauref)
1612
1613	if (tracer.and.(dustbin.ne.0)) then
1614	c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
1615	if (doubleq) then
1616	c call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
1617	c & 'kg.m-2',2,qsurf(1,1))
1618	c call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
1619	c & 'N.m-2',2,qsurf(1,2))
1620	c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
1621	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
1622	c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
1623	c & 'kg.m-2.s-1',2,zdqssed(1,1))
1624	call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
1625	& 'm',3,rdust*ref_r0)
1626	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
1627	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
1628	c call WRITEDIAGFI(ngridmx,'dustN','Dust number',
1629	c & 'part/kg',3,pq(1,1,igcm_dust_number))
1630	#ifdef MESOINI
1631	call WRITEDIAGFI(ngridmx,'dustN','Dust number',
1632	& 'part/kg',3,pq(1,1,igcm_dust_number))
1633	#endif
1634	else
1635	do iq=1,dustbin
1636	write(str2(1:2),'(i2.2)') iq
1637	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
1638	& 'kg/kg',3,zq(1,1,iq))
1639	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
1640	& 'kg.m-2',2,qsurf(1,iq))
1641	end do
1642	endif ! (doubleq)
1643	c if (submicron) then
1644	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
1645	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
1646	c endif ! (submicron)
1647	end if ! (tracer.and.(dustbin.ne.0))
1648
1649	c ----------------------------------------------------------
1650	c ----------------------------------------------------------
1651	c PBL OUTPUS
1652	c ----------------------------------------------------------
1653	c ----------------------------------------------------------
1654
1655
1656	c ----------------------------------------------------------
1657	c Outputs of surface layer
1658	c ----------------------------------------------------------
1659
1660
1661	z_out=1. !0.
1662	if (calltherm .and. (z_out .gt. 0.)) then
1663	call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
1664	& ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
1665
1666	zu2(:)=sqrt(zu(:,1)zu(:,1)+zv(:,1)zv(:,1))
1667	call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
1668	& 'horizontal velocity norm','m/s',
1669	& 2,zu2)
1670
1671	call WRITEDIAGFI(ngridmx,'Teta_out',
1672	& 'potential temperature at z_out','K',
1673	& 2,Teta_out)
1674	call WRITEDIAGFI(ngridmx,'u_out',
1675	& 'horizontal velocity norm at z_out','m/s',
1676	& 2,u_out)
1677	call WRITEDIAGFI(ngridmx,'u*',
1678	& 'friction velocity','m/s',
1679	& 2,ustar)
1680	call WRITEDIAGFI(ngridmx,'teta*',
1681	& 'friction potential temperature','K',
1682	& 2,tstar)
1683	call WRITEDIAGFI(ngrid,'L',
1684	& 'Monin Obukhov length','m',
1685	& 2,L_mo)
1686	else
1687	if((.not. calltherm).and.(z_out .gt. 0.)) then
1688	print*, 'WARNING : no interpolation in surface-layer :'
1689	print*, 'Outputing surface-layer quantities without thermals
1690	& does not make sense'
1691	endif
1692	endif
1693
1694	c ----------------------------------------------------------
1695	c Outputs of thermals
1696	c ----------------------------------------------------------
1697	if (calltherm) then
1698
1699	! call WRITEDIAGFI(ngrid,'dtke',
1700	! & 'tendance tke thermiques','m2/s2',
1701	! & 3,dtke_th)
1702	! call WRITEDIAGFI(ngrid,'d_u_ajs',
1703	! & 'tendance u thermiques','m/s',
1704	! & 3,pdu_th*ptimestep)
1705	! call WRITEDIAGFI(ngrid,'d_v_ajs',
1706	! & 'tendance v thermiques','m/s',
1707	! & 3,pdv_th*ptimestep)
1708	! if (tracer) then
1709	! if (nq .eq. 2) then
1710	! call WRITEDIAGFI(ngrid,'deltaq_th',
1711	! & 'delta q thermiques','kg/kg',
1712	! & 3,ptimestep*pdq_th(:,:,2))
1713	! endif
1714	! endif
1715
1716	call WRITEDIAGFI(ngridmx,'lmax_th',
1717	& 'hauteur du thermique','K',
1718	& 2,lmax_th_out)
1719	call WRITEDIAGFI(ngridmx,'hfmax_th',
1720	& 'maximum TH heat flux','K.m/s',
1721	& 2,hfmax_th)
1722	call WRITEDIAGFI(ngridmx,'wmax_th',
1723	& 'maximum TH vertical velocity','m/s',
1724	& 2,wmax_th)
1725
1726	endif
1727
1728	c ----------------------------------------------------------
1729	c ----------------------------------------------------------
1730	c END OF PBL OUTPUS
1731	c ----------------------------------------------------------
1732	c ----------------------------------------------------------
1733
1734
1735	c ----------------------------------------------------------
1736	c Output in netcdf file "diagsoil.nc" for subterranean
1737	c variables (output every "ecritphy", as for writediagfi)
1738	c ----------------------------------------------------------
1739
1740	! Write soil temperature
1741	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
1742	! & 3,tsoil)
1743	! Write surface temperature
1744	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
1745	! & 2,tsurf)
1746
1747	c ==========================================================
1748	c END OF WRITEDIAGFI
1749	c ==========================================================
1750	#endif
1751
1752	ELSE ! if(ngrid.eq.1)
1753
1754	print,'Ls =',zls180./pi,
1755	& ' tauref(700 Pa) =',tauref
1756	c ----------------------------------------------------------------------
1757	c Output in grads file "g1d" (ONLY when using testphys1d)
1758	c (output at every X physical timestep)
1759	c ----------------------------------------------------------------------
1760	c
1761	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
1762	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
1763	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
1764
1765	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
1766	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
1767	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
1768	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
1769
1770	! THERMALS STUFF 1D
1771
1772	z_out=0.
1773	if (calltherm .and. (z_out .gt. 0.)) then
1774	call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
1775	& ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar)
1776
1777	zu2(:)=sqrt(zu(:,1)zu(:,1)+zv(:,1)zv(:,1))
1778	call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
1779	& 'horizontal velocity norm','m/s',
1780	& 0,zu2)
1781
1782	call WRITEDIAGFI(ngridmx,'Teta_out',
1783	& 'potential temperature at z_out','K',
1784	& 0,Teta_out)
1785	call WRITEDIAGFI(ngridmx,'u_out',
1786	& 'horizontal velocity norm at z_out','m/s',
1787	& 0,u_out)
1788	call WRITEDIAGFI(ngridmx,'u*',
1789	& 'friction velocity','m/s',
1790	& 0,ustar)
1791	call WRITEDIAGFI(ngridmx,'teta*',
1792	& 'friction potential temperature','K',
1793	& 0,tstar)
1794	else
1795	if((.not. calltherm).and.(z_out .gt. 0.)) then
1796	print*, 'WARNING : no interpolation in surface-layer :'
1797	print*, 'Outputing surface-layer quantities without thermals
1798	& does not make sense'
1799	endif
1800	endif
1801
1802	if(calltherm) then
1803
1804	call WRITEDIAGFI(ngridmx,'lmax_th',
1805	& 'hauteur du thermique','point',
1806	& 0,lmax_th_out)
1807	call WRITEDIAGFI(ngridmx,'hfmax_th',
1808	& 'maximum TH heat flux','K.m/s',
1809	& 0,hfmax_th)
1810	call WRITEDIAGFI(ngridmx,'wmax_th',
1811	& 'maximum TH vertical velocity','m/s',
1812	& 0,wmax_th)
1813
1814	co2col(:)=0.
1815	if (tracer) then
1816	do l=1,nlayermx
1817	do ig=1,ngrid
1818	co2col(ig)=co2col(ig)+
1819	& zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
1820	enddo
1821	enddo
1822
1823	end if
1824	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
1825	& ,'kg/m-2',0,co2col)
1826	endif
1827	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
1828	& ,'m/s',1,pw)
1829	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
1830	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
1831	& tsurf)
1832	call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
1833	call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
1834
1835	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
1836	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
1837	! or output in diagfi.nc (for testphys1d)
1838	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
1839	call WRITEDIAGFI(ngridmx,'temp','Temperature',
1840	& 'K',1,zt)
1841
1842	if(tracer) then
1843	c CALL writeg1d(ngrid,1,tau,'tau','SI')
1844	do iq=1,nq
1845	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
1846	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
1847	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
1848	end do
1849	end if
1850
1851	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
1852
1853	do l=2,nlayer-1
1854	tmean=zt(1,l)
1855	if(zt(1,l).ne.zt(1,l-1))
1856	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
1857	zlocal(l)= zlocal(l-1)
1858	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
1859	enddo
1860	zlocal(nlayer)= zlocal(nlayer-1)-
1861	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
1862	& rnew(1,nlayer)*tmean/g
1863
1864	END IF ! if(ngrid.ne.1)
1865
1866	icount=icount+1
1867	RETURN
1868	END

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

Download in other formats:

Original Format