Context Navigation

physiq.F @ 299

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

LMDZ.MARS [associated to previous commit]

21/09/11 == AS

--> Added MESOINI precompiler flag so that all fields needed to initialize the mesoscale model

are being output in the diagfi.nc file. Previously it was made through a separate physiq.F
which needed to be updated every time...

--> This is completely transparent to the casual GCM users and only appears in the WRITEDIAGFI section of physiq.F

File size: 65.9 KB

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

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

Download in other formats:

Original Format