Context Navigation

physiq.F @ 420

Last change on this file since 420 was 420, checked in by tnavarro, 13 years ago

24/11/11 == TN

corrected minor bug in updatereffrad.F : reffdust was not saved

ccn_factor as not correctly used in sedimentation.

It is now initialized in inifis.F, declared in tracer.h and
used in both simpleclouds.F & callsedim.F to update ice radius.

Commented diagfi outputs in aeropacity.F & improvedclouds.F for non scavenging users.

File size: 74.0 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 (micron2/cm3, if photochemistry)
204	REAL surfice(ngridmx,nlayermx) ! ice surface area (micron2/cm3, 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 !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	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	#ifdef MESOSCALE
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
1025	c 7b. Chemical species
1026	c ------------------
1027
1028	#ifndef MESOSCALE
1029	c --------------
1030	c photochemistry :
1031	c --------------
1032	IF (photochem .or. thermochem) then
1033	!NB: Photochemistry includes condensation of H2O2
1034	PRINT*, 'SURFDUST,SURFICE TO BE IMPLEMENTED. YAAAAAARG.'
1035	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
1036	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
1037	$ zdqcloud,zdqscloud,tauref,co2ice,
1038	$ pu,pdu,pv,pdv,surfdust,surfice)
1039
1040	! increment values of tracers:
1041	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1042	! tracers is zero anyways
1043	DO l=1,nlayer
1044	DO ig=1,ngrid
1045	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
1046	ENDDO
1047	ENDDO
1048	ENDDO ! of DO iq=1,nq
1049	! add condensation tendency for H2O2
1050	if (igcm_h2o2.ne.0) then
1051	DO l=1,nlayer
1052	DO ig=1,ngrid
1053	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
1054	& +zdqcloud(ig,l,igcm_h2o2)
1055	ENDDO
1056	ENDDO
1057	endif
1058
1059	! increment surface values of tracers:
1060	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1061	! tracers is zero anyways
1062	DO ig=1,ngrid
1063	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
1064	ENDDO
1065	ENDDO ! of DO iq=1,nq
1066	! add condensation tendency for H2O2
1067	if (igcm_h2o2.ne.0) then
1068	DO ig=1,ngrid
1069	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
1070	& +zdqscloud(ig,igcm_h2o2)
1071	ENDDO
1072	endif
1073
1074	END IF ! of IF (photochem.or.thermochem)
1075	#endif
1076
1077	c 7c. Aerosol particles
1078	c -------------------
1079
1080	c ----------
1081	c Dust devil :
1082	c ----------
1083	IF(callddevil) then
1084	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
1085	& zdqdev,zdqsdev)
1086
1087	if (dustbin.ge.1) then
1088	do iq=1,nq
1089	DO l=1,nlayer
1090	DO ig=1,ngrid
1091	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
1092	ENDDO
1093	ENDDO
1094	enddo
1095	do iq=1,nq
1096	DO ig=1,ngrid
1097	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
1098	ENDDO
1099	enddo
1100	endif ! of if (dustbin.ge.1)
1101
1102	END IF ! of IF (callddevil)
1103
1104	c -------------
1105	c Sedimentation : acts also on water ice
1106	c -------------
1107	IF (sedimentation) THEN
1108	!call zerophys(ngridnlayernq, zdqsed)
1109	zdqsed(1:ngrid,1:nlayer,1:nq)=0
1110	!call zerophys(ngrid*nq, zdqssed)
1111	zdqssed(1:ngrid,1:nq)=0
1112
1113	call callsedim(ngrid,nlayer, ptimestep,
1114	& pplev,zzlev, zzlay, pt, rdust, rice,
1115	& rsedcloud,rhocloud,
1116	& pq, pdq, zdqsed, zdqssed,nq,
1117	& tau,tauscaling)
1118
1119	DO iq=1, nq
1120	DO l=1,nlayer
1121	DO ig=1,ngrid
1122	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1123	ENDDO
1124	ENDDO
1125	ENDDO
1126	DO iq=1, nq
1127	DO ig=1,ngrid
1128	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1129	ENDDO
1130	ENDDO
1131	END IF ! of IF (sedimentation)
1132
1133
1134	c 7d. Updates
1135	c ---------
1136
1137	DO iq=1, nq
1138	DO ig=1,ngrid
1139
1140	c ---------------------------------
1141	c Updating tracer budget on surface
1142	c ---------------------------------
1143	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1144
1145	ENDDO ! (ig)
1146	ENDDO ! (iq)
1147
1148	endif ! of if (tracer)
1149
1150	#ifndef MESOSCALE
1151	c-----------------------------------------------------------------------
1152	c 8. THERMOSPHERE CALCULATION
1153	c-----------------------------------------------------------------------
1154
1155	if (callthermos) then
1156	call thermosphere(pplev,pplay,dist_sol,
1157	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1158	& pt,pq,pu,pv,pdt,pdq,
1159	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1160
1161	DO l=1,nlayer
1162	DO ig=1,ngrid
1163	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1164	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1165	& +zdteuv(ig,l)
1166	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1167	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1168	DO iq=1, nq
1169	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1170	ENDDO
1171	ENDDO
1172	ENDDO
1173
1174	endif ! of if (callthermos)
1175	#endif
1176
1177	c-----------------------------------------------------------------------
1178	c 9. Surface and sub-surface soil temperature
1179	c-----------------------------------------------------------------------
1180	c
1181	c
1182	c 9.1 Increment Surface temperature:
1183	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1184
1185	DO ig=1,ngrid
1186	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1187	ENDDO
1188
1189	c Prescribe a cold trap at south pole (except at high obliquity !!)
1190	c Temperature at the surface is set there to be the temperature
1191	c corresponding to equilibrium temperature between phases of CO2
1192
1193
1194	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
1195	#ifndef MESOSCALE
1196	if (caps.and.(obliquit.lt.27.)) then
1197	! NB: Updated surface pressure, at grid point 'ngrid', is
1198	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1199	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1200	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1201	endif
1202	#endif
1203	c -------------------------------------------------------------
1204	c Change of surface albedo in case of ground frost
1205	c everywhere except on the north permanent cap and in regions
1206	c covered by dry ice.
1207	c ALWAYS PLACE these lines after newcondens !!!
1208	c -------------------------------------------------------------
1209	do ig=1,ngrid
1210	if ((co2ice(ig).eq.0).and.
1211	& (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
1212	albedo(ig,1) = albedo_h2o_ice
1213	albedo(ig,2) = albedo_h2o_ice
1214	c write(,) "frost thickness", qsurf(ig,igcm_h2o_ice)
1215	c write(,) "physiq.F frost :"
1216	c & ,lati(ig)180./pi, long(ig)180./pi
1217	endif
1218	enddo ! of do ig=1,ngrid
1219	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
1220
1221	c
1222	c 9.2 Compute soil temperatures and subsurface heat flux:
1223	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1224	IF (callsoil) THEN
1225	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1226	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
1227	ENDIF
1228
1229
1230	c-----------------------------------------------------------------------
1231	c 10. Write output files
1232	c ----------------------
1233
1234	c -------------------------------
1235	c Dynamical fields incrementation
1236	c -------------------------------
1237	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1238	! temperature, zonal and meridional wind
1239	DO l=1,nlayer
1240	DO ig=1,ngrid
1241	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1242	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1243	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1244	ENDDO
1245	ENDDO
1246
1247	! tracers
1248	DO iq=1, nq
1249	DO l=1,nlayer
1250	DO ig=1,ngrid
1251	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1252	ENDDO
1253	ENDDO
1254	ENDDO
1255
1256	! surface pressure
1257	DO ig=1,ngrid
1258	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1259	ENDDO
1260
1261	! pressure
1262	DO l=1,nlayer
1263	DO ig=1,ngrid
1264	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1265	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1266	ENDDO
1267	ENDDO
1268
1269	! Density
1270	DO l=1,nlayer
1271	DO ig=1,ngrid
1272	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1273	ENDDO
1274	ENDDO
1275
1276	! Potential Temperature
1277
1278	DO ig=1,ngridmx
1279	DO l=1,nlayermx
1280	zh(ig,l) = zt(ig,l)(zplay(ig,l)/zplev(ig,1))*rcp
1281	ENDDO
1282	ENDDO
1283
1284
1285	c Compute surface stress : (NB: z0 is a common in surfdat.h)
1286	c DO ig=1,ngrid
1287	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
1288	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1289	c ENDDO
1290
1291	c Sum of fluxes in solar spectral bands (for output only)
1292	DO ig=1,ngrid
1293	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1294	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1295	ENDDO
1296	c ***** TEST ****************************************************
1297	ztim1 = 999
1298	DO l=1,nlayer
1299	DO ig=1,ngrid
1300	if (pt(ig,l).lt.ztim1) then
1301	ztim1 = pt(ig,l)
1302	igmin = ig
1303	lmin = l
1304	end if
1305	ENDDO
1306	ENDDO
1307	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1308	write(,) 'PHYSIQ: stability WARNING :'
1309	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1310	& 'ig l =', igmin, lmin
1311	end if
1312	c *******************************************************************
1313
1314	c ---------------------
1315	c Outputs to the screen
1316	c ---------------------
1317
1318	IF (lwrite) THEN
1319	PRINT*,'Global diagnostics for the physics'
1320	PRINT,'Variables and their increments x and dx/dt dt'
1321	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1322	WRITE(*,'(2f10.5,2f15.5)')
1323	s tsurf(igout),zdtsurf(igout)*ptimestep,
1324	s pplev(igout,1),pdpsrf(igout)*ptimestep
1325	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1326	WRITE(*,'(i4,6f10.5)') (l,
1327	s pu(igout,l),pdu(igout,l)*ptimestep,
1328	s pv(igout,l),pdv(igout,l)*ptimestep,
1329	s pt(igout,l),pdt(igout,l)*ptimestep,
1330	s l=1,nlayer)
1331	ENDIF ! of IF (lwrite)
1332
1333	IF (ngrid.NE.1) THEN
1334
1335	#ifndef MESOSCALE
1336	c -------------------------------------------------------------------
1337	c Writing NetCDF file "RESTARTFI" at the end of the run
1338	c -------------------------------------------------------------------
1339	c Note: 'restartfi' is stored just before dynamics are stored
1340	c in 'restart'. Between now and the writting of 'restart',
1341	c there will have been the itau=itau+1 instruction and
1342	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1343	c thus we store for time=time+dtvr
1344
1345	IF(lastcall) THEN
1346	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1347	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1348	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1349	. ptimestep,pday,
1350	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
1351	. area,albedodat,inertiedat,zmea,zstd,zsig,
1352	. zgam,zthe)
1353	ENDIF
1354	#endif
1355
1356	c -------------------------------------------------------------------
1357	c Calculation of diagnostic variables written in both stats and
1358	c diagfi files
1359	c -------------------------------------------------------------------
1360
1361	if (tracer) then
1362	if (water) then
1363
1364	mtot(:)=0
1365	icetot(:)=0
1366	rave(:)=0
1367	tauTES(:)=0
1368	do ig=1,ngrid
1369	do l=1,nlayermx
1370	mtot(ig) = mtot(ig) +
1371	& zq(ig,l,igcm_h2o_vap) *
1372	& (pplev(ig,l) - pplev(ig,l+1)) / g
1373	icetot(ig) = icetot(ig) +
1374	& zq(ig,l,igcm_h2o_ice) *
1375	& (pplev(ig,l) - pplev(ig,l+1)) / g
1376	rave(ig) = rave(ig) +
1377	& zq(ig,l,igcm_h2o_ice) *
1378	& (pplev(ig,l) - pplev(ig,l+1)) / g *
1379	& rice(ig,l) * (1.+nuice_ref)
1380	c Computing abs optical depth at 825 cm-1 in each
1381	c layer to simulate NEW TES retrieval
1382	Qabsice = min(
1383	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
1384	& )
1385	opTES(ig,l)= 0.75 * Qabsice *
1386	& zq(ig,l,igcm_h2o_ice) *
1387	& (pplev(ig,l) - pplev(ig,l+1)) / g
1388	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1389	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1390	enddo
1391	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
1392	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1393	enddo
1394
1395	endif ! of if (water)
1396	endif ! of if (tracer)
1397
1398	c -----------------------------------------------------------------
1399	c WSTATS: Saving statistics
1400	c -----------------------------------------------------------------
1401	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1402	c which can later be used to make the statistic files of the run:
1403	c "stats") only possible in 3D runs !
1404
1405	IF (callstats) THEN
1406
1407	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1408	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1409	call wstats(ngrid,"co2ice","CO2 ice cover",
1410	& "kg.m-2",2,co2ice)
1411	call wstats(ngrid,"fluxsurf_lw",
1412	& "Thermal IR radiative flux to surface","W.m-2",2,
1413	& fluxsurf_lw)
1414	call wstats(ngrid,"fluxsurf_sw",
1415	& "Solar radiative flux to surface","W.m-2",2,
1416	& fluxsurf_sw_tot)
1417	call wstats(ngrid,"fluxtop_lw",
1418	& "Thermal IR radiative flux to space","W.m-2",2,
1419	& fluxtop_lw)
1420	call wstats(ngrid,"fluxtop_sw",
1421	& "Solar radiative flux to space","W.m-2",2,
1422	& fluxtop_sw_tot)
1423	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1424	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1425	call wstats(ngrid,"v","Meridional (North-South) wind",
1426	& "m.s-1",3,zv)
1427	call wstats(ngrid,"w","Vertical (down-up) wind",
1428	& "m.s-1",3,pw)
1429	call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
1430	call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
1431	c call wstats(ngrid,"q2",
1432	c & "Boundary layer eddy kinetic energy",
1433	c & "m2.s-2",3,q2)
1434	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1435	c & emis)
1436	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
1437	c & 2,zstress)
1438	c call wstats(ngrid,"sw_htrt","sw heat.rate",
1439	c & "W.m-2",3,zdtsw)
1440	c call wstats(ngrid,"lw_htrt","lw heat.rate",
1441	c & "W.m-2",3,zdtlw)
1442
1443	if (tracer) then
1444	if (water) then
1445	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1446	c & *mugaz/mmol(igcm_h2o_vap)
1447	c call wstats(ngrid,"vmr_h2ovapor",
1448	c & "H2O vapor volume mixing ratio","mol/mol",
1449	c & 3,vmr)
1450	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1451	c & *mugaz/mmol(igcm_h2o_ice)
1452	c call wstats(ngrid,"vmr_h2oice",
1453	c & "H2O ice volume mixing ratio","mol/mol",
1454	c & 3,vmr)
1455	call wstats(ngrid,"h2o_ice_s",
1456	& "surface h2o_ice","kg/m2",
1457	& 2,qsurf(1,igcm_h2o_ice))
1458	c call wstats(ngrid,'albedo',
1459	c & 'albedo',
1460	c & '',2,albedo(1:ngridmx,1))
1461	call wstats(ngrid,"mtot",
1462	& "total mass of water vapor","kg/m2",
1463	& 2,mtot)
1464	call wstats(ngrid,"icetot",
1465	& "total mass of water ice","kg/m2",
1466	& 2,icetot)
1467	c call wstats(ngrid,"reffice",
1468	c & "Mean reff","m",
1469	c & 2,rave)
1470	c call wstats(ngrid,"rice",
1471	c & "Ice particle size","m",
1472	c & 3,rice)
1473	c If activice is true, tauTES is computed in aeropacity.F;
1474	if (.not.activice) then
1475	call wstats(ngrid,"tauTESap",
1476	& "tau abs 825 cm-1","",
1477	& 2,tauTES)
1478	endif
1479
1480	endif ! of if (water)
1481
1482	if (thermochem.or.photochem) then
1483	do iq=1,nq
1484	if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
1485	. (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
1486	. (noms(iq).eq."h2").or.
1487	. (noms(iq).eq."o3")) then
1488	do l=1,nlayer
1489	do ig=1,ngrid
1490	vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
1491	end do
1492	end do
1493	call wstats(ngrid,"vmr_"//trim(noms(iq)),
1494	. "Volume mixing ratio","mol/mol",3,vmr)
1495	endif
1496	! do ig = 1,ngrid
1497	! colden(ig,iq) = 0. !FL
1498	! end do
1499	! do l=1,nlayer !FL
1500	! do ig=1,ngrid !FL
1501	! colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq) !FL
1502	! $ *(pplev(ig,l)-pplev(ig,l+1)) !FL
1503	! $ 6.022e22/(mmol(iq)g) !FL
1504	! end do !FL
1505	! end do !FL
1506	! call wstats(ngrid,"c_"//trim(noms(iq)), !FL
1507	! $ "column","mol cm-2",2,colden(1,iq)) !FL
1508	end do
1509	end if ! of if (thermochem.or.photochem)
1510
1511	end if ! of if (tracer)
1512
1513	IF(lastcall) THEN
1514	write (,) "Writing stats..."
1515	call mkstats(ierr)
1516	ENDIF
1517
1518	ENDIF !if callstats
1519
1520	c (Store EOF for Mars Climate database software)
1521	IF (calleofdump) THEN
1522	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1523	ENDIF
1524
1525
1526	#ifdef MESOSCALE
1527	!!!
1528	!!! OUTPUT FIELDS
1529	!!!
1530	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
1531	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
1532	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
1533	!! JF
1534	TAU_lay(:)=tau(:,1)!!true opacity (not a reference like tauref)
1535	IF (igcm_dust_mass .ne. 0) THEN
1536	qsurfice_dust(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
1537	ENDIF
1538	IF (igcm_h2o_ice .ne. 0) THEN
1539	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
1540	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1541	. * mugaz / mmol(igcm_h2o_ice)
1542	ENDIF
1543	!! Dust quantity integration along the vertical axe
1544	dustot(:)=0
1545	do ig=1,ngrid
1546	do l=1,nlayermx
1547	dustot(ig) = dustot(ig) +
1548	& zq(ig,l,igcm_dust_mass)
1549	& * (pplev(ig,l) - pplev(ig,l+1)) / g
1550	enddo
1551	enddo
1552	c AUTOMATICALLY GENERATED FROM REGISTRY
1553	#include "fill_save.inc"
1554	#else
1555
1556	c ==========================================================
1557	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
1558	c any variable for diagnostic (output with period
1559	c "ecritphy", set in "run.def")
1560	c ==========================================================
1561	c WRITEDIAGFI can ALSO be called from any other subroutines
1562	c for any variables !!
1563	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1564	c & emis)
1565	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
1566	! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
1567	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1568	& tsurf)
1569	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1570	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
1571	& co2ice)
1572
1573	c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
1574	c & "K",2,zt(1,7))
1575	c call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1576	c & fluxsurf_lw)
1577	c call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1578	c & fluxsurf_sw_tot)
1579	c call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1580	c & fluxtop_lw)
1581	c call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1582	c & fluxtop_sw_tot)
1583	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1584	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1585	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1586	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1587	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1588	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1589	! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1590	c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
1591	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
1592	c & zstress)
1593	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
1594	c & 'w.m-2',3,zdtsw)
1595	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
1596	c & 'w.m-2',3,zdtlw)
1597	c CALL WRITEDIAGFI(ngridmx,'tauTESap',
1598	c & 'tau abs 825 cm-1',
1599	c & '',2,tauTES)
1600
1601	#ifdef MESOINI
1602	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1603	& emis)
1604	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1605	call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
1606	& "K",3,tsoil)
1607	call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
1608	& "K",3,inertiedat)
1609	#endif
1610
1611
1612	c ----------------------------------------------------------
1613	c Outputs of the CO2 cycle
1614	c ----------------------------------------------------------
1615
1616	if (tracer.and.(igcm_co2.ne.0)) then
1617	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
1618	! & "kg/kg",2,zq(1,1,igcm_co2))
1619	! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
1620	! & "kg/kg",3,zq(1,1,igcm_co2))
1621
1622	! Compute co2 column
1623	co2col(:)=0
1624	do l=1,nlayermx
1625	do ig=1,ngrid
1626	co2col(ig)=co2col(ig)+
1627	& zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
1628	enddo
1629	enddo
1630	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
1631	& co2col)
1632	!!!!! FL
1633	! do iq = 1,nq
1634	! if (noms(iq) .ne. "dust_mass" .and.
1635	! $ noms(iq) .ne. "dust_number") then
1636	! call writediagfi(ngrid,"c_"//trim(noms(iq)),
1637	! $ "column","mol cm-2",2,colden(1,iq))
1638	! end if
1639	! end do
1640	!!!!! FL
1641	endif ! of if (tracer.and.(igcm_co2.ne.0))
1642
1643	c ----------------------------------------------------------
1644	c Outputs of the water cycle
1645	c ----------------------------------------------------------
1646	if (tracer) then
1647	if (water) then
1648
1649	#ifdef MESOINI
1650	!!!! waterice = q01, voir readmeteo.F90
1651	call WRITEDIAGFI(ngridmx,'q01',noms(igcm_h2o_ice),
1652	& 'kg/kg',3,
1653	& zq(1:ngridmx,1:nlayermx,igcm_h2o_ice))
1654	!!!! watervapor = q02, voir readmeteo.F90
1655	call WRITEDIAGFI(ngridmx,'q02',noms(igcm_h2o_vap),
1656	& 'kg/kg',3,
1657	& zq(1:ngridmx,1:nlayermx,igcm_h2o_vap))
1658	!!!! surface waterice qsurf02 (voir readmeteo)
1659	call WRITEDIAGFI(ngridmx,'qsurf02','surface tracer',
1660	& 'kg.m-2',2,
1661	& qsurf(1:ngridmx,igcm_h2o_ice))
1662	#endif
1663
1664	CALL WRITEDIAGFI(ngridmx,'mtot',
1665	& 'total mass of water vapor',
1666	& 'kg/m2',2,mtot)
1667	CALL WRITEDIAGFI(ngridmx,'icetot',
1668	& 'total mass of water ice',
1669	& 'kg/m2',2,icetot)
1670	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1671	c & *mugaz/mmol(igcm_h2o_ice)
1672	c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
1673	c & 'mol/mol',3,vmr)
1674	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1675	c & *mugaz/mmol(igcm_h2o_vap)
1676	c call WRITEDIAGFI(ngridmx,'vmr_h2ovap','h2o vap vmr',
1677	c & 'mol/mol',3,vmr)
1678	c CALL WRITEDIAGFI(ngridmx,'reffice',
1679	c & 'Mean reff',
1680	c & 'm',2,rave)
1681	c call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
1682	c & 'm',3,rice)
1683	c If activice is true, tauTES is computed in aeropacity.F;
1684	if (.not.activice) then
1685	CALL WRITEDIAGFI(ngridmx,'tauTESap',
1686	& 'tau abs 825 cm-1',
1687	& '',2,tauTES)
1688	endif
1689	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
1690	& 'surface h2o_ice',
1691	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
1692
1693	if (caps) then
1694	do ig=1,ngridmx
1695	if (watercaptag(ig)) watercapflag(ig) = 1
1696	enddo
1697	c CALL WRITEDIAGFI(ngridmx,'watercaptag',
1698	c & 'Ice water caps',
1699	c & '',2,watercapflag)
1700	endif
1701	c CALL WRITEDIAGFI(ngridmx,'albedo',
1702	c & 'albedo',
1703	c & '',2,albedo(1:ngridmx,1))
1704	endif !(water)
1705
1706
1707	if (water.and..not.photochem) then
1708	iq=nq
1709	c write(str2(1:2),'(i2.2)') iq
1710	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
1711	c & 'kg.m-2',2,zdqscloud(1,iq))
1712	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
1713	c & 'kg/kg',3,zdqchim(1,1,iq))
1714	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
1715	c & 'kg/kg',3,zdqdif(1,1,iq))
1716	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
1717	c & 'kg/kg',3,zdqadj(1,1,iq))
1718	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
1719	c & 'kg/kg',3,zdqc(1,1,iq))
1720	endif !(water.and..not.photochem)
1721	endif
1722
1723	c ----------------------------------------------------------
1724	c Outputs of the dust cycle
1725	c ----------------------------------------------------------
1726
1727	c call WRITEDIAGFI(ngridmx,'tauref',
1728	c & 'Dust ref opt depth','NU',2,tauref)
1729
1730	if (tracer.and.(dustbin.ne.0)) then
1731	c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
1732	if (doubleq) then
1733	c call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
1734	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
1735	c call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
1736	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
1737	c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
1738	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
1739	c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
1740	c & 'kg.m-2.s-1',2,zdqssed(1,1))
1741	c call WRITEDIAGFI(ngridmx,'dqsdif','diffusion',
1742	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
1743	c call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
1744	c & 'm',3,rdust*ref_r0)
1745	c call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
1746	c & 'kg/kg',3,pq(1,1,igcm_dust_mass))
1747	c call WRITEDIAGFI(ngridmx,'dustN','Dust number',
1748	c & 'part/kg',3,pq(1,1,igcm_dust_number))
1749	#ifdef MESOINI
1750	call WRITEDIAGFI(ngridmx,'dustN','Dust number',
1751	& 'part/kg',3,pq(1,1,igcm_dust_number))
1752	#endif
1753	else
1754	do iq=1,dustbin
1755	write(str2(1:2),'(i2.2)') iq
1756	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
1757	& 'kg/kg',3,zq(1,1,iq))
1758	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
1759	& 'kg.m-2',2,qsurf(1,iq))
1760	end do
1761	endif ! (doubleq)
1762
1763	if (scavenging) then
1764	call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr',
1765	& 'kg/kg',3,pq(1,1,igcm_ccn_mass))
1766	call WRITEDIAGFI(ngridmx,'ccnN','CCN number',
1767	& 'part/kg',3,pq(1,1,igcm_ccn_number))
1768	endif ! (scavenging)
1769
1770	c if (submicron) then
1771	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
1772	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
1773	c endif ! (submicron)
1774	end if ! (tracer.and.(dustbin.ne.0))
1775
1776	c ----------------------------------------------------------
1777	c ----------------------------------------------------------
1778	c PBL OUTPUS
1779	c ----------------------------------------------------------
1780	c ----------------------------------------------------------
1781
1782
1783	c ----------------------------------------------------------
1784	c Outputs of surface layer
1785	c ----------------------------------------------------------
1786
1787
1788	z_out=0.
1789	if (calltherm .and. (z_out .gt. 0.)) then
1790	call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
1791	& ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
1792
1793	zu2(:)=sqrt(zu(:,1)zu(:,1)+zv(:,1)zv(:,1))
1794	call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
1795	& 'horizontal velocity norm','m/s',
1796	& 2,zu2)
1797
1798	call WRITEDIAGFI(ngridmx,'Teta_out',
1799	& 'potential temperature at z_out','K',
1800	& 2,Teta_out)
1801	call WRITEDIAGFI(ngridmx,'u_out',
1802	& 'horizontal velocity norm at z_out','m/s',
1803	& 2,u_out)
1804	call WRITEDIAGFI(ngridmx,'u*',
1805	& 'friction velocity','m/s',
1806	& 2,ustar)
1807	call WRITEDIAGFI(ngridmx,'teta*',
1808	& 'friction potential temperature','K',
1809	& 2,tstar)
1810	call WRITEDIAGFI(ngrid,'L',
1811	& 'Monin Obukhov length','m',
1812	& 2,L_mo)
1813	else
1814	if((.not. calltherm).and.(z_out .gt. 0.)) then
1815	print*, 'WARNING : no interpolation in surface-layer :'
1816	print*, 'Outputing surface-layer quantities without thermals
1817	& does not make sense'
1818	endif
1819	endif
1820
1821	c ----------------------------------------------------------
1822	c Outputs of thermals
1823	c ----------------------------------------------------------
1824	if (calltherm) then
1825
1826	! call WRITEDIAGFI(ngrid,'dtke',
1827	! & 'tendance tke thermiques','m2/s2',
1828	! & 3,dtke_th)
1829	! call WRITEDIAGFI(ngrid,'d_u_ajs',
1830	! & 'tendance u thermiques','m/s',
1831	! & 3,pdu_th*ptimestep)
1832	! call WRITEDIAGFI(ngrid,'d_v_ajs',
1833	! & 'tendance v thermiques','m/s',
1834	! & 3,pdv_th*ptimestep)
1835	! if (tracer) then
1836	! if (nq .eq. 2) then
1837	! call WRITEDIAGFI(ngrid,'deltaq_th',
1838	! & 'delta q thermiques','kg/kg',
1839	! & 3,ptimestep*pdq_th(:,:,2))
1840	! endif
1841	! endif
1842
1843	call WRITEDIAGFI(ngridmx,'lmax_th',
1844	& 'hauteur du thermique','K',
1845	& 2,lmax_th_out)
1846	call WRITEDIAGFI(ngridmx,'hfmax_th',
1847	& 'maximum TH heat flux','K.m/s',
1848	& 2,hfmax_th)
1849	call WRITEDIAGFI(ngridmx,'wmax_th',
1850	& 'maximum TH vertical velocity','m/s',
1851	& 2,wmax_th)
1852
1853	endif
1854
1855	c ----------------------------------------------------------
1856	c ----------------------------------------------------------
1857	c END OF PBL OUTPUS
1858	c ----------------------------------------------------------
1859	c ----------------------------------------------------------
1860
1861
1862	c ----------------------------------------------------------
1863	c Output in netcdf file "diagsoil.nc" for subterranean
1864	c variables (output every "ecritphy", as for writediagfi)
1865	c ----------------------------------------------------------
1866
1867	! Write soil temperature
1868	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
1869	! & 3,tsoil)
1870	! Write surface temperature
1871	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
1872	! & 2,tsurf)
1873
1874	c ==========================================================
1875	c END OF WRITEDIAGFI
1876	c ==========================================================
1877	#endif
1878
1879	ELSE ! if(ngrid.eq.1)
1880
1881	print,'Ls =',zls180./pi,
1882	& ' tauref(700 Pa) =',tauref
1883	c ----------------------------------------------------------------------
1884	c Output in grads file "g1d" (ONLY when using testphys1d)
1885	c (output at every X physical timestep)
1886	c ----------------------------------------------------------------------
1887	c
1888	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
1889	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
1890	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
1891
1892	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
1893	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
1894	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
1895	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
1896
1897	! THERMALS STUFF 1D
1898
1899	z_out=0.
1900	if (calltherm .and. (z_out .gt. 0.)) then
1901	call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
1902	& ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
1903
1904	zu2(:)=sqrt(zu(:,1)zu(:,1)+zv(:,1)zv(:,1))
1905	call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
1906	& 'horizontal velocity norm','m/s',
1907	& 0,zu2)
1908
1909	call WRITEDIAGFI(ngridmx,'Teta_out',
1910	& 'potential temperature at z_out','K',
1911	& 0,Teta_out)
1912	call WRITEDIAGFI(ngridmx,'u_out',
1913	& 'horizontal velocity norm at z_out','m/s',
1914	& 0,u_out)
1915	call WRITEDIAGFI(ngridmx,'u*',
1916	& 'friction velocity','m/s',
1917	& 0,ustar)
1918	call WRITEDIAGFI(ngridmx,'teta*',
1919	& 'friction potential temperature','K',
1920	& 0,tstar)
1921	else
1922	if((.not. calltherm).and.(z_out .gt. 0.)) then
1923	print*, 'WARNING : no interpolation in surface-layer :'
1924	print*, 'Outputing surface-layer quantities without thermals
1925	& does not make sense'
1926	endif
1927	endif
1928
1929	if(calltherm) then
1930
1931	call WRITEDIAGFI(ngridmx,'lmax_th',
1932	& 'hauteur du thermique','point',
1933	& 0,lmax_th_out)
1934	call WRITEDIAGFI(ngridmx,'hfmax_th',
1935	& 'maximum TH heat flux','K.m/s',
1936	& 0,hfmax_th)
1937	call WRITEDIAGFI(ngridmx,'wmax_th',
1938	& 'maximum TH vertical velocity','m/s',
1939	& 0,wmax_th)
1940
1941	co2col(:)=0.
1942	if (tracer) then
1943	do l=1,nlayermx
1944	do ig=1,ngrid
1945	co2col(ig)=co2col(ig)+
1946	& zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
1947	enddo
1948	enddo
1949
1950	end if
1951	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
1952	& ,'kg/m-2',0,co2col)
1953	endif
1954	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
1955	& ,'m/s',1,pw)
1956	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
1957	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
1958	& tsurf)
1959	call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
1960	call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
1961
1962	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
1963	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
1964	call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
1965	! or output in diagfi.nc (for testphys1d)
1966	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
1967	call WRITEDIAGFI(ngridmx,'temp','Temperature',
1968	& 'K',1,zt)
1969
1970	if(tracer) then
1971	c CALL writeg1d(ngrid,1,tau,'tau','SI')
1972	do iq=1,nq
1973	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
1974	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
1975	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
1976	end do
1977	if (doubleq) then
1978	call WRITEDIAGFI(ngridmx,'rdust','rdust',
1979	& 'm',1,rdust)
1980	endif
1981	end if
1982
1983	cccccccccccccccccccccc scavenging outputs 1D TN ccccccccccccccccccc
1984	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1985	IF (scavenging) THEN
1986	CALL WRITEDIAGFI(ngridmx,'tauTESap',
1987	& 'tau abs 825 cm-1',
1988	& '',1,tauTES)
1989
1990	h2o_tot = qsurf(1,igcm_h2o_ice)
1991	do l=1,nlayer
1992	h2o_tot = h2o_tot +
1993	& (zq(1,l,igcm_h2o_vap) + zq(1,l,igcm_h2o_ice))
1994	& * (pplev(1,l) - pplev(1,l+1)) / g
1995	end do
1996
1997
1998	do iq=1,nq
1999	call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_s',
2000	& trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
2001	end do
2002
2003
2004	call WRITEDIAGFI(ngridmx,'dqssed q','sedimentation q',
2005	& 'kg.m-2.s-1',0,zdqssed(:,igcm_dust_number))
2006	call WRITEDIAGFI(ngridmx,'dqsdif q','diffusion q',
2007	& 'kg.m-2.s-1',0,zdqsdif(:,igcm_dust_number))
2008
2009	call WRITEDIAGFI(ngridmx,'dqssed N','sedimentation N',
2010	& 'kg.m-2.s-1',0,zdqssed(:,igcm_dust_number))
2011	call WRITEDIAGFI(ngridmx,'dqsdif N','diffusion N',
2012	& 'kg.m-2.s-1',0,zdqsdif(:,igcm_dust_number))
2013
2014	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
2015	& rice)
2016	ENDIF
2017
2018	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2019	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2020
2021
2022	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
2023
2024	do l=2,nlayer-1
2025	tmean=zt(1,l)
2026	if(zt(1,l).ne.zt(1,l-1))
2027	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
2028	zlocal(l)= zlocal(l-1)
2029	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
2030	enddo
2031	zlocal(nlayer)= zlocal(nlayer-1)-
2032	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
2033	& rnew(1,nlayer)*tmean/g
2034
2035	END IF ! if(ngrid.ne.1)
2036
2037	icount=icount+1
2038	RETURN
2039	END

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

Download in other formats:

Original Format