Context Navigation

physiq.F @ 350

Last change on this file since 350 was 334, checked in by aslmd, 14 years ago

LMDZ.MARS:

28/10/11 == FL + AS

ADDED THE NEW FRAMEWORK FOR PHOTOCHEMISTRY
This is not the last version. A new rewritten version of calchim.F [using LAPACK] is yet to be committed.
--> A new version for photochemistry routines : *_B.F no longer exists, new routines in aeronomars
D 333 libf/aeronomars/photochemist_B.F
D 333 libf/aeronomars/init_chimie_B.F
A 0 libf/aeronomars/read_phototable.F
M 333 libf/aeronomars/calchim.F
A 0 libf/aeronomars/photochemistry.F
M 333 libf/aeronomars/chimiedata.h
--> Changed calls to calchim and watercloud [surfdust and surfice needed] in physiq.F
--> Left commented code for outputs in physiq.F [search for FL]
--> Added settings which works for 35 levels in inidissip.F according to FL. Commented for the moment.
--> Checked compilation and run, looks fine. Note that 'jmars.20101220' is needed.

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

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

Download in other formats:

Original Format