Context Navigation

physiq.F @ 485

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

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

Download in other formats:

Original Format