Context Navigation

physiq.F @ 267

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

--- AC 03/08/2011 ---
M 265 libf/phymars/physiq.F
<> Added a PBL section for outputs, with a call to SL outputs via surflayer_interpol and thermals outputs

A 0 libf/phymars/surflayer_interpol.F
<> New subroutine to interpolate horizontal velocity norm and potential temperature in the surface layer.

THIS ROUTINE IS NOT VALIDATED YET. IT IS TURNED OFF BY DEFAULT AND IS HERE FOR DEVELOPMENT PURPOSES ONLY FOR NOW.

M 265 libf/phymars/vdif_cd.F
<> Important modification to the Reynolds formula : due to a confusion of symbols used in the litterature, a wrong simplification of

numerical values had been implicitely done... It is now corrected and yields better results for the mixed layer temperature.

M 265 libf/phymars/vdifc.F
<> Cosmetic modif (added comment for more clarity)

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

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

Download in other formats:

Original Format