Context Navigation

physiq.F @ 286

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

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

Download in other formats:

Original Format