Context Navigation

physiq.F @ 1246

Last change on this file since 1246 was 1246, checked in by aslmd, 11 years ago
LMDZ.MARS. Made number of scatterers a free dimension not in need to be prescribe at compiling time. Instead it must be set in callphys.def. See README for further information about this commit.
File size: 95.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
10	use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2,
11	& igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
12	& igcm_ccn_mass, igcm_ccn_number,
13	& igcm_dust_mass, igcm_dust_number, igcm_h2o2,
14	& nuice_ref, rho_ice, rho_dust, ref_r0
15	use comsoil_h, only: inertiedat, ! soil thermal inertia
16	& tsoil, nsoilmx ! number of subsurface layers
17	use comgeomfi_h, only: long, lati, area,
18	& sinlon, coslon, sinlat, coslat
19	use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
20	& zthe, z0, albedo_h2o_ice,
21	& frost_albedo_threshold,
22	& tsurf, co2ice, emis,
23	& capcal, fluxgrd, qsurf
24	use comsaison_h, only: dist_sol, declin, mu0, fract
25	use slope_mod, only: theta_sl, psi_sl
26	use conc_mod, only: rnew, cpnew, mmean
27	use control_mod, only: iphysiq, day_step, ecritstart
28	use dimradmars_mod, only: tauscaling, aerosol,
29	& dtrad, fluxrad_sky, fluxrad, albedo,
30	& naerkind
31	use turb_mod, only: q2, wstar, ustar, sensibFlux,
32	& zmax_th, hfmax_th, turb_resolved
33	use planete_h, only: aphelie, periheli, year_day, peri_day,
34	& obliquit
35	USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad, daysec
36	#ifdef MESOSCALE
37	use comsoil_h, only: mlayer,layer
38	use surfdat_h, only: z0_default
39	use comm_wrf
40	#else
41	use phyredem, only: physdem0, physdem1
42	use eofdump_mod, only: eofdump
43	#endif
44
45	IMPLICIT NONE
46	c=======================================================================
47	c
48	c subject:
49	c --------
50	c
51	c Organisation of the physical parametrisations of the LMD
52	c martian atmospheric general circulation model.
53	c
54	c The GCM can be run without or with tracer transport
55	c depending on the value of Logical "tracer" in file "callphys.def"
56	c Tracers may be water vapor, ice OR chemical species OR dust particles
57	c
58	c SEE comments in initracer.F about numbering of tracer species...
59	c
60	c It includes:
61	c
62	c 1. Initialization:
63	c 1.1 First call initializations
64	c 1.2 Initialization for every call to physiq
65	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
66	c 2. Compute radiative transfer tendencies
67	c (longwave and shortwave) for CO2 and aerosols.
68	c 3. Gravity wave and subgrid scale topography drag :
69	c 4. Vertical diffusion (turbulent mixing):
70	c 5. Convective adjustment
71	c 6. Condensation and sublimation of carbon dioxide.
72	c 7. TRACERS :
73	c 7a. water and water ice
74	c 7b. call for photochemistry when tracers are chemical species
75	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
76	c 7d. updates (CO2 pressure variations, surface budget)
77	c 8. Contribution to tendencies due to thermosphere
78	c 9. Surface and sub-surface temperature calculations
79	c 10. Write outputs :
80	c - "startfi", "histfi" (if it's time)
81	c - Saving statistics (if "callstats = .true.")
82	c - Dumping eof (if "calleofdump = .true.")
83	c - Output any needed variables in "diagfi"
84	c 11. Diagnostic: mass conservation of tracers
85	c
86	c author:
87	c -------
88	c Frederic Hourdin 15/10/93
89	c Francois Forget 1994
90	c Christophe Hourdin 02/1997
91	c Subroutine completly rewritten by F.Forget (01/2000)
92	c Introduction of the photochemical module: S. Lebonnois (11/2002)
93	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
94	c Water ice clouds: Franck Montmessin (update 06/2003)
95	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
96	c Nb: See callradite.F for more information.
97	c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
98	c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
99	c
100	c arguments:
101	c ----------
102	c
103	c input:
104	c ------
105	c ecri period (in dynamical timestep) to write output
106	c ngrid Size of the horizontal grid.
107	c All internal loops are performed on that grid.
108	c nlayer Number of vertical layers.
109	c nq Number of advected fields
110	c firstcall True at the first call
111	c lastcall True at the last call
112	c pday Number of days counted from the North. Spring
113	c equinoxe.
114	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
115	c ptimestep timestep (s)
116	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
117	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
118	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
119	c pu(ngrid,nlayer) u component of the wind (ms-1)
120	c pv(ngrid,nlayer) v component of the wind (ms-1)
121	c pt(ngrid,nlayer) Temperature (K)
122	c pq(ngrid,nlayer,nq) Advected fields
123	c pudyn(ngrid,nlayer) \|
124	c pvdyn(ngrid,nlayer) \| Dynamical temporal derivative for the
125	c ptdyn(ngrid,nlayer) \| corresponding variables
126	c pqdyn(ngrid,nlayer,nq) \|
127	c pw(ngrid,?) vertical velocity
128	c
129	c output:
130	c -------
131	c
132	c pdu(ngrid,nlayer) \|
133	c pdv(ngrid,nlayer) \| Temporal derivative of the corresponding
134	c pdt(ngrid,nlayer) \| variables due to physical processes.
135	c pdq(ngrid,nlayer,nq) \|
136	c pdpsrf(ngrid) \|
137	c tracerdyn call tracer in dynamical part of GCM ?
138
139	c
140	c=======================================================================
141	c
142	c 0. Declarations :
143	c ------------------
144
145	#include "dimensions.h"
146	#include "dimphys.h"
147	!#include "comgeomfi.h"
148	!#include "surfdat.h"
149	!#include "comsoil.h"
150	!#include "comdiurn.h"
151	#include "callkeys.h"
152	!#include "comsaison.h"
153	!#include "control.h"
154	!#include "dimradmars.h"
155	#include "comg1d.h"
156	!#include "tracer.h"
157	#include "nlteparams.h"
158	#include "comvert.h"
159
160	#include "chimiedata.h"
161	#include "param.h"
162	#include "param_v4.h"
163	!#include "conc.h"
164
165	#include "netcdf.inc"
166
167	!#include "slope.h"
168
169	c Arguments :
170	c -----------
171
172	c inputs:
173	c -------
174	INTEGER ngrid,nlayer,nq
175	REAL ptimestep
176	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
177	REAL pphi(ngrid,nlayer)
178	REAL pu(ngrid,nlayer),pv(ngrid,nlayer)
179	REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
180	REAL pw(ngrid,nlayer) !Mars pvervel transmit par dyn3d
181	REAL zh(ngrid,nlayer) ! potential temperature (K)
182	LOGICAL firstcall,lastcall
183
184	REAL pday
185	REAL ptime
186	logical tracerdyn
187
188	c outputs:
189	c --------
190	c physical tendencies
191	REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer)
192	REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
193	REAL pdpsrf(ngrid) ! surface pressure tendency
194
195
196	c Local saved variables:
197	c ----------------------
198	INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
199	INTEGER,SAVE :: icount ! counter of calls to physiq during the run.
200
201	c Variables used by the water ice microphysical scheme:
202	REAL rice(ngrid,nlayer) ! Water ice geometric mean radius (m)
203	REAL nuice(ngrid,nlayer) ! Estimated effective variance
204	! of the size distribution
205	real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius
206	real rhocloud(ngrid,nlayer) ! Cloud density (kg.m-3)
207	REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
208	REAL surfice(ngrid,nlayer) ! ice surface area (m2/m3, if photochemistry)
209	REAL inertiesoil(ngrid,nsoilmx) ! Time varying subsurface
210	! thermal inertia (J.s-1/2.m-2.K-1)
211	! (used only when tifeedback=.true.)
212
213	c Variables used by the slope model
214	REAL sl_ls, sl_lct, sl_lat
215	REAL sl_tau, sl_alb, sl_the, sl_psi
216	REAL sl_fl0, sl_flu
217	REAL sl_ra, sl_di0
218	REAL sky
219
220	REAL,PARAMETER :: stephan = 5.67e-08 ! Stephan Boltzman constant
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,tapphys
230
231	REAL fluxsurf_lw(ngrid) !incident LW (IR) surface flux (W.m-2)
232	REAL fluxsurf_sw(ngrid,2) !incident SW (solar) surface flux (W.m-2)
233	REAL fluxtop_lw(ngrid) !Outgoing LW (IR) flux to space (W.m-2)
234	REAL fluxtop_sw(ngrid,2) !Outgoing SW (solar) flux to space (W.m-2)
235	REAL tauref(ngrid) ! Reference column optical depth at odpref
236	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
237	REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
238	! AS: TBD: this one should be in a module !
239	REAL zls ! solar longitude (rad)
240	REAL zday ! date (time since Ls=0, in martian days)
241	REAL zzlay(ngrid,nlayer) ! altitude at the middle of the layers
242	REAL zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
243	! REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
244
245	c Tendancies due to various processes:
246	REAL dqsurf(ngrid,nq)
247	REAL zdtlw(ngrid,nlayer) ! (K/s)
248	REAL zdtsw(ngrid,nlayer) ! (K/s)
249	! REAL cldtlw(ngrid,nlayer) ! (K/s) LW heating rate for clear area
250	! REAL cldtsw(ngrid,nlayer) ! (K/s) SW heating rate for clear area
251	REAL zdtnirco2(ngrid,nlayer) ! (K/s)
252	REAL zdtnlte(ngrid,nlayer) ! (K/s)
253	REAL zdtsurf(ngrid) ! (K/s)
254	REAL zdtcloud(ngrid,nlayer)
255	REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! (m.s-2)
256	REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid) ! (K/s)
257	REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! (m.s-2)
258	REAL zdhadj(ngrid,nlayer) ! (K/s)
259	REAL zdtgw(ngrid,nlayer) ! (K/s)
260	REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! (m.s-2)
261	REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid)
262	REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
263
264	REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
265	REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
266	REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
267	REAL zdqadj(ngrid,nlayer,nq)
268	REAL zdqc(ngrid,nlayer,nq)
269	REAL zdqcloud(ngrid,nlayer,nq)
270	REAL zdqscloud(ngrid,nq)
271	REAL zdqchim(ngrid,nlayer,nq)
272	REAL zdqschim(ngrid,nq)
273
274	REAL zdteuv(ngrid,nlayer) ! (K/s)
275	REAL zdtconduc(ngrid,nlayer) ! (K/s)
276	REAL zdumolvis(ngrid,nlayer)
277	REAL zdvmolvis(ngrid,nlayer)
278	real zdqmoldiff(ngrid,nlayer,nq)
279
280	c Local variable for local intermediate calcul:
281	REAL zflubid(ngrid)
282	REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
283	REAL zdum1(ngrid,nlayer)
284	REAL zdum2(ngrid,nlayer)
285	REAL ztim1,ztim2,ztim3, z1,z2
286	REAL ztime_fin
287	REAL zdh(ngrid,nlayer)
288	INTEGER length
289	PARAMETER (length=100)
290
291	c local variables only used for diagnostic (output in file "diagfi" or "stats")
292	c -----------------------------------------------------------------------------
293	REAL ps(ngrid), zt(ngrid,nlayer)
294	REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
295	REAL zq(ngrid,nlayer,nq)
296	REAL fluxtop_sw_tot(ngrid), fluxsurf_sw_tot(ngrid)
297	character*2 str2
298	! character*5 str5
299	real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
300	real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
301	integer igmin, lmin
302	logical tdiag
303
304	real co2col(ngrid) ! CO2 column
305	! pplev and pplay are dynamical inputs and must not be modified in the physics.
306	! instead, use zplay and zplev :
307	REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
308	! REAL zstress(ngrid),cd
309	real tmean, zlocal(nlayer)
310	real rho(ngrid,nlayer) ! density
311	real vmr(ngrid,nlayer) ! volume mixing ratio
312	real rhopart(ngrid,nlayer) ! number density of a given species
313	real colden(ngrid,nq) ! vertical column of tracers
314	REAL mtot(ngrid) ! Total mass of water vapor (kg/m2)
315	REAL icetot(ngrid) ! Total mass of water ice (kg/m2)
316	REAL Nccntot(ngrid) ! Total number of ccn (nbr/m2)
317	REAL Mccntot(ngrid) ! Total mass of ccn (kg/m2)
318	REAL rave(ngrid) ! Mean water ice effective radius (m)
319	REAL opTES(ngrid,nlayer) ! abs optical depth at 825 cm-1
320	REAL tauTES(ngrid) ! column optical depth at 825 cm-1
321	REAL Qabsice ! Water ice absorption coefficient
322	REAL taucloudtes(ngrid) ! Cloud opacity at infrared
323	! reference wavelength using
324	! Qabs instead of Qext
325	! (direct comparison with TES)
326
327	REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
328	REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
329	REAL ndust(ngrid,nlayer) ! true n dust (kg/kg)
330	REAL qdust(ngrid,nlayer) ! true q dust (kg/kg)
331	REAL nccn(ngrid,nlayer) ! true n ccn (kg/kg)
332	REAL qccn(ngrid,nlayer) ! true q ccn (kg/kg)
333
334	c Test 1d/3d scavenging
335	real h2otot(ngrid)
336	REAL satu(ngrid,nlayer) ! satu ratio for output
337	REAL zqsat(ngrid,nlayer) ! saturation
338
339	REAL,SAVE :: time_phys
340
341	! Added for new NLTE scheme
342
343	real co2vmr_gcm(ngrid,nlayer)
344	real n2vmr_gcm(ngrid,nlayer)
345	real ovmr_gcm(ngrid,nlayer)
346	real covmr_gcm(ngrid,nlayer)
347	integer ierr_nlte
348	real*8 varerr
349
350	c Variables for PBL
351	REAL zz1(ngrid)
352	REAL lmax_th_out(ngrid)
353	REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
354	REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
355	INTEGER lmax_th(ngrid),dimout,n_out,n
356	CHARACTER(50) zstring
357	REAL dtke_th(ngrid,nlayer+1)
358	REAL zcdv(ngrid), zcdh(ngrid)
359	REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
360	REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
361	REAL u_out1(ngrid)
362	REAL T_out1(ngrid)
363	REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters]
364	REAL tstar(ngrid) ! friction velocity and friction potential temp
365	REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid)
366	! REAL zu2(ngrid)
367
368	c=======================================================================
369
370	c 1. Initialisation:
371	c -----------------
372
373	c 1.1 Initialisation only at first call
374	c ---------------------------------------
375	IF (firstcall) THEN
376
377	c variables set to 0
378	c ~~~~~~~~~~~~~~~~~~
379	aerosol(:,:,:)=0
380	dtrad(:,:)=0
381
382	#ifndef MESOSCALE
383	fluxrad(:)=0
384	wstar(:)=0.
385	#endif
386
387	c read startfi
388	c ~~~~~~~~~~~~
389	#ifndef MESOSCALE
390	! GCM. Read netcdf initial physical parameters.
391	CALL phyetat0 ("startfi.nc",0,0,
392	& nsoilmx,ngrid,nlayer,nq,
393	& day_ini,time_phys,
394	& tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling)
395
396	if (pday.ne.day_ini) then
397	write(,) "PHYSIQ: ERROR: bad synchronization between ",
398	& "physics and dynamics"
399	write(,) "dynamics day: ",pday
400	write(,) "physics day: ",day_ini
401	stop
402	endif
403
404	write (,) 'In physiq day_ini =', day_ini
405
406	#else
407	! MESOSCALE. Supposedly everything is already set in modules.
408	! So we just check. And we calculate day_ini + two param in dyn3d/control_mod.
409	print*,"check: rad,cpp,g,r,rcp,daysec"
410	print*,rad,cpp,g,r,rcp,daysec
411	PRINT*,'check: tsurf ',tsurf(1),tsurf(ngrid)
412	PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngrid,nsoilmx)
413	PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
414	PRINT*,'check: midlayer,layer ', mlayer(:),layer(:)
415	PRINT*,'check: tracernames ', noms
416	PRINT*,'check: emis ',emis(1),emis(ngrid)
417	PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayermx+1)
418	PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngrid,nq)
419	PRINT*,'check: co2 ',co2ice(1),co2ice(ngrid)
420	day_step=daysec/ptimestep
421	PRINT*,'Call to LMD physics:',day_step,' per Martian day'
422	iphysiq=ptimestep
423	day_ini = pday
424	#endif
425
426	c initialize tracers
427	c ~~~~~~~~~~~~~~~~~~
428	tracerdyn=tracer
429	IF (tracer) THEN
430	CALL initracer(ngrid,nq,qsurf)
431	ENDIF ! end tracer
432
433	c Initialize albedo and orbital calculation
434	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
435	CALL surfini(ngrid,co2ice,qsurf,albedo)
436	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
437
438	c initialize soil
439	c ~~~~~~~~~~~~~~~
440	IF (callsoil) THEN
441	c Thermal inertia feedback:
442	IF (tifeedback) THEN
443	CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
444	CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
445	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
446	ELSE
447	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
448	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
449	ENDIF ! of IF (tifeedback)
450	ELSE
451	PRINT*,
452	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
453	DO ig=1,ngrid
454	capcal(ig)=1.e5
455	fluxgrd(ig)=0.
456	ENDDO
457	ENDIF
458	icount=1
459
460	#ifndef MESOSCALE
461	c Initialize thermospheric parameters
462	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
463
464	if (callthermos) then
465	if(solvarmod.eq.0) call param_read
466	if(solvarmod.eq.1) call param_read_e107
467	endif
468	#endif
469	c Initialize R and Cp as constant
470
471	if (.not.callthermos .and. .not.photochem) then
472	do l=1,nlayer
473	do ig=1,ngrid
474	rnew(ig,l)=r
475	cpnew(ig,l)=cpp
476	mmean(ig,l)=mugaz
477	enddo
478	enddo
479	endif
480
481	if(callnlte.and.nltemodel.eq.2) call nlte_setup
482	if(callnirco2.and.nircorr.eq.1) call NIR_leedat
483	if(thermochem) call chemthermos_readini
484
485	IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
486	write(,)"physiq: water_param Surface water ice albedo:",
487	. albedo_h2o_ice
488	ENDIF
489
490	#ifndef MESOSCALE
491	if (callslope) call getslopes(ngrid,phisfi)
492
493	if (ngrid.ne.1) then ! no need to create a restart file in 1d
494	call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq,
495	. ptimestep,pday,time_phys,area,
496	. albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
497	endif
498	#endif
499
500	ENDIF ! (end of "if firstcall")
501
502
503	c ---------------------------------------------------
504	c 1.2 Initializations done at every physical timestep:
505	c ---------------------------------------------------
506	c
507
508	c Initialize various variables
509	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
510	pdv(:,:)=0
511	pdu(:,:)=0
512	pdt(:,:)=0
513	pdq(:,:,:)=0
514	pdpsrf(:)=0
515	zflubid(:)=0
516	zdtsurf(:)=0
517	dqsurf(:,:)=0
518	igout=ngrid/2+1
519
520
521	zday=pday+ptime ! compute time, in sols (and fraction thereof)
522
523	c Compute Solar Longitude (Ls) :
524	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
525	if (season) then
526	call solarlong(zday,zls)
527	else
528	call solarlong(float(day_ini),zls)
529	end if
530
531	c Initialize pressure levels
532	c ~~~~~~~~~~~~~~~~~~
533	zplev(:,:) = pplev(:,:)
534	zplay(:,:) = pplay(:,:)
535	ps(:) = pplev(:,1)
536
537
538	c Compute geopotential at interlayers
539	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
540	c ponderation des altitudes au niveau des couches en dp/p
541
542	DO l=1,nlayer
543	DO ig=1,ngrid
544	zzlay(ig,l)=pphi(ig,l)/g
545	ENDDO
546	ENDDO
547	DO ig=1,ngrid
548	zzlev(ig,1)=0.
549	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
550	ENDDO
551	DO l=2,nlayer
552	DO ig=1,ngrid
553	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
554	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
555	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
556	ENDDO
557	ENDDO
558
559
560	! Potential temperature calculation not the same in physiq and dynamic
561
562	c Compute potential temperature
563	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
564	DO l=1,nlayer
565	DO ig=1,ngrid
566	zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
567	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
568	ENDDO
569	ENDDO
570
571	#ifndef MESOSCALE
572	c-----------------------------------------------------------------------
573	c 1.2.5 Compute mean mass, cp, and R
574	c --------------------------------
575
576	if(photochem.or.callthermos) then
577	call concentrations(ngrid,nlayer,nq,
578	& zplay,pt,pdt,pq,pdq,ptimestep)
579	endif
580	#endif
581	c-----------------------------------------------------------------------
582	c 2. Compute radiative tendencies :
583	c------------------------------------
584
585
586	IF (callrad) THEN
587	IF( MOD(icount-1,iradia).EQ.0) THEN
588
589	c Local Solar zenith angle
590	c ~~~~~~~~~~~~~~~~~~~~~~~~
591	CALL orbite(zls,dist_sol,declin)
592
593	IF(diurnal) THEN
594	ztim1=SIN(declin)
595	ztim2=COS(declin)COS(2.pi*(zday-.5))
596	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
597
598	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
599	s ztim1,ztim2,ztim3, mu0,fract)
600
601	ELSE
602	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
603	ENDIF
604
605	c NLTE cooling from CO2 emission
606	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
607	IF(callnlte) then
608	if(nltemodel.eq.0.or.nltemodel.eq.1) then
609	CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
610	else if(nltemodel.eq.2) then
611	co2vmr_gcm(1:ngrid,1:nlayer)=
612	& pq(1:ngrid,1:nlayer,igcm_co2)*
613	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
614	n2vmr_gcm(1:ngrid,1:nlayer)=
615	& pq(1:ngrid,1:nlayer,igcm_n2)*
616	& mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
617	covmr_gcm(1:ngrid,1:nlayer)=
618	& pq(1:ngrid,1:nlayer,igcm_co)*
619	& mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
620	ovmr_gcm(1:ngrid,1:nlayer)=
621	& pq(1:ngrid,1:nlayer,igcm_o)*
622	& mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
623
624	CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
625	$ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
626	$ ovmr_gcm, zdtnlte,ierr_nlte,varerr )
627	if(ierr_nlte.gt.0) then
628	write(,)
629	$ 'WARNING: nlte_tcool output with error message',
630	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
631	write(,)'I will continue anyway'
632	endif
633
634	zdtnlte(1:ngrid,1:nlayer)=
635	& zdtnlte(1:ngrid,1:nlayer)/86400.
636	endif
637	else
638	zdtnlte(:,:)=0.
639	endif
640
641	c Find number of layers for LTE radiation calculations
642	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
643	& CALL nlthermeq(ngrid,nlayer,zplev,zplay)
644
645	c Note: Dustopacity.F has been transferred to callradite.F
646
647	c Call main radiative transfer scheme
648	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
649	c Transfer through CO2 (except NIR CO2 absorption)
650	c and aerosols (dust and water ice)
651
652	c Radiative transfer
653	c ------------------
654	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
655	$ emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
656	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
657	$ tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice,
658	$ nuice,co2ice)
659
660	c Outputs for basic check (middle of domain)
661	c ------------------------------------------
662	write(*,'("Ls =",f11.6," check lat =",f10.6,
663	& " lon =",f11.6)')
664	& zls180./pi,lati(igout)180/pi,long(igout)*180/pi
665	write(*,'(" tauref(",f4.0," Pa) =",f9.6,
666	& " tau(",f4.0," Pa) =",f9.6)')
667	& odpref,tauref(igout),
668	& odpref,tau(igout,1)*odpref/zplev(igout,1)
669	c ---------------------------------------------------------
670	c Call slope parameterization for direct and scattered flux
671	c ---------------------------------------------------------
672	IF(callslope) THEN
673	print *, 'Slope scheme is on and computing...'
674	DO ig=1,ngrid
675	sl_the = theta_sl(ig)
676	IF (sl_the .ne. 0.) THEN
677	ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
678	DO l=1,2
679	sl_lct = ptime24. + 180.long(ig)/pi/15.
680	sl_ra = pi*(1.0-sl_lct/12.)
681	sl_lat = 180.*lati(ig)/pi
682	sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
683	sl_alb = albedo(ig,l)
684	sl_psi = psi_sl(ig)
685	sl_fl0 = fluxsurf_sw(ig,l)
686	sl_di0 = 0.
687	if (mu0(ig) .gt. 0.) then
688	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
689	sl_di0 = sl_di0*1370./dist_sol/dist_sol
690	sl_di0 = sl_di0/ztim1
691	sl_di0 = fluxsurf_sw(ig,l)*sl_di0
692	endif
693	! you never know (roundup concern...)
694	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
695	!!!!!!!!!!!!!!!!!!!!!!!!!!
696	CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
697	& sl_tau, sl_alb, sl_the, sl_psi,
698	& sl_di0, sl_fl0, sl_flu )
699	!!!!!!!!!!!!!!!!!!!!!!!!!!
700	fluxsurf_sw(ig,l) = sl_flu
701	ENDDO
702	!!! compute correction on IR flux as well
703	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
704	fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
705	ENDIF
706	ENDDO
707	ENDIF
708
709	c CO2 near infrared absorption
710	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
711	zdtnirco2(:,:)=0
712	if (callnirco2) then
713	call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
714	. mu0,fract,declin, zdtnirco2)
715	endif
716
717	c Radiative flux from the sky absorbed by the surface (W.m-2)
718	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
719	DO ig=1,ngrid
720	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
721	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
722	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
723	ENDDO
724
725
726	c Net atmospheric radiative heating rate (K.s-1)
727	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
728	IF(callnlte) THEN
729	CALL blendrad(ngrid, nlayer, zplay,
730	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
731	ELSE
732	DO l=1,nlayer
733	DO ig=1,ngrid
734	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
735	& +zdtnirco2(ig,l)
736	ENDDO
737	ENDDO
738	ENDIF
739
740	ENDIF ! of if(mod(icount-1,iradia).eq.0)
741
742	c Transformation of the radiative tendencies:
743	c -------------------------------------------
744
745	c Net radiative surface flux (W.m-2)
746	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
747	c
748	DO ig=1,ngrid
749	zplanck(ig)=tsurf(ig)*tsurf(ig)
750	zplanck(ig)=emis(ig)*
751	$ stephanzplanck(ig)zplanck(ig)
752	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
753	IF(callslope) THEN
754	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
755	fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
756	ENDIF
757	ENDDO
758
759	DO l=1,nlayer
760	DO ig=1,ngrid
761	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
762	ENDDO
763	ENDDO
764
765	ENDIF ! of IF (callrad)
766
767	c-----------------------------------------------------------------------
768	c 3. Gravity wave and subgrid scale topography drag :
769	c -------------------------------------------------
770
771
772	IF(calllott)THEN
773
774	CALL calldrag_noro(ngrid,nlayer,ptimestep,
775	& zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
776
777	DO l=1,nlayer
778	DO ig=1,ngrid
779	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
780	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
781	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
782	ENDDO
783	ENDDO
784	ENDIF
785
786	c-----------------------------------------------------------------------
787	c 4. Vertical diffusion (turbulent mixing):
788	c -----------------------------------------
789
790	IF (calldifv) THEN
791
792	DO ig=1,ngrid
793	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
794	ENDDO
795
796	zdum1(:,:)=0
797	zdum2(:,:)=0
798	do l=1,nlayer
799	do ig=1,ngrid
800	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
801	enddo
802	enddo
803
804	c ----------------------
805	c Treatment of a special case : using new surface layer (Richardson based)
806	c without using the thermals in gcm and mesoscale can yield problems in
807	c weakly unstable situations when winds are near to 0. For those cases, we add
808	c a unit subgrid gustiness. Remember that thermals should be used we using the
809	c Richardson based surface layer model.
810	IF ( .not.calltherm
811	. .and. callrichsl
812	. .and. .not.turb_resolved) THEN
813	DO ig=1, ngrid
814	IF (zh(ig,1) .lt. tsurf(ig)) THEN
815	wstar(ig)=1.
816	hfmax_th(ig)=0.2
817	ELSE
818	wstar(ig)=0.
819	hfmax_th(ig)=0.
820	ENDIF
821	ENDDO
822	ENDIF
823	c ----------------------
824
825	IF (tke_heat_flux .ne. 0.) THEN
826	zz1(:)=(pt(:,1)+pdt(:,1)ptimestep)(r/g)*
827	& (-alog(zplay(:,1)/zplev(:,1)))
828	pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
829	ENDIF
830
831	c Calling vdif (Martian version WITH CO2 condensation)
832	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
833	$ ptimestep,capcal,lwrite,
834	$ zplay,zplev,zzlay,zzlev,z0,
835	$ pu,pv,zh,pq,tsurf,emis,qsurf,
836	$ zdum1,zdum2,zdh,pdq,zflubid,
837	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
838	& zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux)
839
840
841	DO ig=1,ngrid
842	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
843	ENDDO
844
845	IF (.not.turb_resolved) THEN
846	DO l=1,nlayer
847	DO ig=1,ngrid
848	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
849	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
850	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
851
852	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
853	ENDDO
854	ENDDO
855
856	if (tracer) then
857	DO iq=1, nq
858	DO l=1,nlayer
859	DO ig=1,ngrid
860	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
861	ENDDO
862	ENDDO
863	ENDDO
864	DO iq=1, nq
865	DO ig=1,ngrid
866	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
867	ENDDO
868	ENDDO
869	end if ! of if (tracer)
870	ELSE
871	write (,) '******************************************'
872	write (,) '** LES mode: the difv part is only used to'
873	write (,) '** - provide HFX and UST to the dynamics'
874	write (,) '** - update TSURF'
875	write (,) '******************************************'
876	!! Specific treatment for lifting in turbulent-resolving mode (AC)
877	IF (lifting .and. doubleq) THEN
878	!! lifted dust is injected in the first layer.
879	!! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF.
880	!! => lifted dust is not incremented before the sedimentation step.
881	zdqdif(1:ngrid,1,1:nq)=0.
882	zdqdif(1:ngrid,1,igcm_dust_number) =
883	. -zdqsdif(1:ngrid,igcm_dust_number)
884	zdqdif(1:ngrid,1,igcm_dust_mass) =
885	. -zdqsdif(1:ngrid,igcm_dust_mass)
886	zdqdif(1:ngrid,2:nlayermx,1:nq) = 0.
887	DO iq=1, nq
888	IF ((iq .ne. igcm_dust_mass)
889	& .and. (iq .ne. igcm_dust_number)) THEN
890	zdqsdif(:,iq)=0.
891	ENDIF
892	ENDDO
893	ELSE
894	zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
895	zdqsdif(1:ngrid,1:nq) = 0.
896	ENDIF
897	ENDIF
898	ELSE
899	DO ig=1,ngrid
900	zdtsurf(ig)=zdtsurf(ig)+
901	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
902	ENDDO
903	IF (turb_resolved) THEN
904	write(,) 'Turbulent-resolving mode !'
905	write(,) 'Please set calldifv to T in callphys.def'
906	STOP
907	ENDIF
908	ENDIF ! of IF (calldifv)
909
910	c-----------------------------------------------------------------------
911	c 5. Thermals :
912	c -----------------------------
913
914	if(calltherm .and. .not.turb_resolved) then
915
916	call calltherm_interface(ngrid,nlayer,nq,
917	$ tracer,igcm_co2,
918	$ zzlev,zzlay,
919	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
920	$ zplay,zplev,pphi,zpopsk,
921	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
922	$ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
923
924	DO l=1,nlayer
925	DO ig=1,ngrid
926	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
927	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
928	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
929	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
930	ENDDO
931	ENDDO
932
933	DO ig=1,ngrid
934	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
935	ENDDO
936
937	if (tracer) then
938	DO iq=1,nq
939	DO l=1,nlayer
940	DO ig=1,ngrid
941	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
942	ENDDO
943	ENDDO
944	ENDDO
945	endif
946
947	lmax_th_out(:)=real(lmax_th(:))
948
949	else !of if calltherm
950	lmax_th(:)=0
951	wstar(:)=0.
952	hfmax_th(:)=0.
953	lmax_th_out(:)=0.
954	end if
955
956	c-----------------------------------------------------------------------
957	c 5. Dry convective adjustment:
958	c -----------------------------
959
960	IF(calladj) THEN
961
962	DO l=1,nlayer
963	DO ig=1,ngrid
964	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
965	ENDDO
966	ENDDO
967	zduadj(:,:)=0
968	zdvadj(:,:)=0
969	zdhadj(:,:)=0
970
971	CALL convadj(ngrid,nlayer,nq,ptimestep,
972	$ zplay,zplev,zpopsk,lmax_th,
973	$ pu,pv,zh,pq,
974	$ pdu,pdv,zdh,pdq,
975	$ zduadj,zdvadj,zdhadj,
976	$ zdqadj)
977
978
979	DO l=1,nlayer
980	DO ig=1,ngrid
981	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
982	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
983	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
984
985	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
986	ENDDO
987	ENDDO
988
989	if(tracer) then
990	DO iq=1, nq
991	DO l=1,nlayer
992	DO ig=1,ngrid
993	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
994	ENDDO
995	ENDDO
996	ENDDO
997	end if
998	ENDIF ! of IF(calladj)
999
1000
1001
1002	c-----------------------------------------------------------------------
1003	c 6. Specific parameterizations for tracers
1004	c: -----------------------------------------
1005
1006	if (tracer) then
1007
1008	c 6a. Water and ice
1009	c ---------------
1010
1011	c ---------------------------------------
1012	c Water ice condensation in the atmosphere
1013	c ----------------------------------------
1014	IF (water) THEN
1015
1016	call watercloud(ngrid,nlayer,ptimestep,
1017	& zplev,zplay,pdpsrf,zzlay, pt,pdt,
1018	& pq,pdq,zdqcloud,zdtcloud,
1019	& nq,tau,tauscaling,rdust,rice,nuice,
1020	& rsedcloud,rhocloud)
1021
1022	c Temperature variation due to latent heat release
1023	if (activice) then
1024	pdt(1:ngrid,1:nlayer) =
1025	& pdt(1:ngrid,1:nlayer) +
1026	& zdtcloud(1:ngrid,1:nlayer)
1027	endif
1028
1029	! increment water vapour and ice atmospheric tracers tendencies
1030	pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
1031	& pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
1032	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
1033	pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1034	& pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1035	& zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
1036
1037	! increment dust and ccn masses and numbers
1038	! We need to check that we have Nccn & Ndust > 0
1039	! This is due to single precision rounding problems
1040	if (microphys) then
1041	pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1042	& pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1043	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
1044	pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1045	& pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1046	& zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
1047	where (pq(:,:,igcm_ccn_mass) +
1048	& ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1049	pdq(:,:,igcm_ccn_mass) =
1050	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1051	pdq(:,:,igcm_ccn_number) =
1052	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1053	end where
1054	where (pq(:,:,igcm_ccn_number) +
1055	& ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1056	pdq(:,:,igcm_ccn_mass) =
1057	& - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1058	pdq(:,:,igcm_ccn_number) =
1059	& - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1060	end where
1061	endif
1062
1063	if (scavenging) then
1064	pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
1065	& pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
1066	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
1067	pdq(1:ngrid,1:nlayer,igcm_dust_number) =
1068	& pdq(1:ngrid,1:nlayer,igcm_dust_number) +
1069	& zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
1070	where (pq(:,:,igcm_dust_mass) +
1071	& ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1072	pdq(:,:,igcm_dust_mass) =
1073	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1074	pdq(:,:,igcm_dust_number) =
1075	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1076	end where
1077	where (pq(:,:,igcm_dust_number) +
1078	& ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1079	pdq(:,:,igcm_dust_mass) =
1080	& - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1081	pdq(:,:,igcm_dust_number) =
1082	& - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1083	end where
1084	endif ! of if scavenging
1085
1086
1087	END IF ! of IF (water)
1088
1089	c 6b. Aerosol particles
1090	c -------------------
1091
1092	c ----------
1093	c Dust devil :
1094	c ----------
1095	IF(callddevil) then
1096	call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
1097	& zdqdev,zdqsdev)
1098
1099	if (dustbin.ge.1) then
1100	do iq=1,nq
1101	DO l=1,nlayer
1102	DO ig=1,ngrid
1103	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
1104	ENDDO
1105	ENDDO
1106	enddo
1107	do iq=1,nq
1108	DO ig=1,ngrid
1109	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
1110	ENDDO
1111	enddo
1112	endif ! of if (dustbin.ge.1)
1113
1114	END IF ! of IF (callddevil)
1115
1116	c -------------
1117	c Sedimentation : acts also on water ice
1118	c -------------
1119	IF (sedimentation) THEN
1120	zdqsed(1:ngrid,1:nlayer,1:nq)=0
1121	zdqssed(1:ngrid,1:nq)=0
1122
1123	call callsedim(ngrid,nlayer, ptimestep,
1124	& zplev,zzlev, zzlay, pt, pdt, rdust, rice,
1125	& rsedcloud,rhocloud,
1126	& pq, pdq, zdqsed, zdqssed,nq,
1127	& tau,tauscaling)
1128
1129	DO iq=1, nq
1130	DO l=1,nlayer
1131	DO ig=1,ngrid
1132	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1133	ENDDO
1134	ENDDO
1135	ENDDO
1136	DO iq=1, nq
1137	DO ig=1,ngrid
1138	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1139	ENDDO
1140	ENDDO
1141	END IF ! of IF (sedimentation)
1142
1143	c Add lifted dust to tendancies after sedimentation in the LES (AC)
1144	IF (turb_resolved) THEN
1145	DO iq=1, nq
1146	DO l=1,nlayer
1147	DO ig=1,ngrid
1148	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
1149	ENDDO
1150	ENDDO
1151	ENDDO
1152	DO iq=1, nq
1153	DO ig=1,ngrid
1154	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
1155	ENDDO
1156	ENDDO
1157	ENDIF
1158
1159	c
1160	c 6c. Chemical species
1161	c ------------------
1162
1163	#ifndef MESOSCALE
1164	c --------------
1165	c photochemistry :
1166	c --------------
1167	IF (photochem .or. thermochem) then
1168
1169	! dust and ice surface area
1170	call surfacearea(ngrid, nlayer, naerkind,
1171	$ ptimestep, zplay, zzlay,
1172	$ pt, pq, pdq, nq,
1173	$ rdust, rice, tau, tauscaling,
1174	$ surfdust, surfice)
1175	! call photochemistry
1176	call calchim(ngrid,nlayer,nq,
1177	& ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
1178	$ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
1179	$ zdqcloud,zdqscloud,tauref,co2ice,
1180	$ pu,pdu,pv,pdv,surfdust,surfice)
1181
1182	! increment values of tracers:
1183	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1184	! tracers is zero anyways
1185	DO l=1,nlayer
1186	DO ig=1,ngrid
1187	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
1188	ENDDO
1189	ENDDO
1190	ENDDO ! of DO iq=1,nq
1191
1192	! add condensation tendency for H2O2
1193	if (igcm_h2o2.ne.0) then
1194	DO l=1,nlayer
1195	DO ig=1,ngrid
1196	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
1197	& +zdqcloud(ig,l,igcm_h2o2)
1198	ENDDO
1199	ENDDO
1200	endif
1201
1202	! increment surface values of tracers:
1203	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1204	! tracers is zero anyways
1205	DO ig=1,ngrid
1206	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
1207	ENDDO
1208	ENDDO ! of DO iq=1,nq
1209
1210	! add condensation tendency for H2O2
1211	if (igcm_h2o2.ne.0) then
1212	DO ig=1,ngrid
1213	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
1214	& +zdqscloud(ig,igcm_h2o2)
1215	ENDDO
1216	endif
1217
1218	END IF ! of IF (photochem.or.thermochem)
1219	#endif
1220
1221	c 6d. Updates
1222	c ---------
1223
1224	DO iq=1, nq
1225	DO ig=1,ngrid
1226
1227	c ---------------------------------
1228	c Updating tracer budget on surface
1229	c ---------------------------------
1230	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1231
1232	ENDDO ! (ig)
1233	ENDDO ! (iq)
1234
1235	endif ! of if (tracer)
1236
1237	#ifndef MESOSCALE
1238	c-----------------------------------------------------------------------
1239	c 7. THERMOSPHERE CALCULATION
1240	c-----------------------------------------------------------------------
1241
1242	if (callthermos) then
1243	call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
1244	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1245	& pt,pq,pu,pv,pdt,pdq,
1246	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1247
1248	DO l=1,nlayer
1249	DO ig=1,ngrid
1250	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1251	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1252	& +zdteuv(ig,l)
1253	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1254	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1255	DO iq=1, nq
1256	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1257	ENDDO
1258	ENDDO
1259	ENDDO
1260
1261	endif ! of if (callthermos)
1262	#endif
1263
1264	c-----------------------------------------------------------------------
1265	c 8. Carbon dioxide condensation-sublimation:
1266	c (should be the last atmospherical physical process to be computed)
1267	c -------------------------------------------
1268
1269	IF (tituscap) THEN
1270	!!! get the actual co2 seasonal cap from Titus observations
1271	CALL geticecover( ngrid, 180.*zls/pi,
1272	. 180.long/pi, 180.lati/pi, co2ice )
1273	co2ice = co2ice * 10000.
1274	ENDIF
1275
1276
1277	pdpsrf(:) = 0
1278
1279	IF (callcond) THEN
1280	CALL newcondens(ngrid,nlayer,nq,ptimestep,
1281	$ capcal,zplay,zplev,tsurf,pt,
1282	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
1283	$ co2ice,albedo,emis,
1284	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
1285	$ fluxsurf_sw,zls)
1286
1287	DO l=1,nlayer
1288	DO ig=1,ngrid
1289	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
1290	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
1291	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
1292	ENDDO
1293	ENDDO
1294	DO ig=1,ngrid
1295	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
1296	ENDDO
1297
1298	IF (tracer) THEN
1299	DO iq=1, nq
1300	DO l=1,nlayer
1301	DO ig=1,ngrid
1302	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1303	ENDDO
1304	ENDDO
1305	ENDDO
1306	ENDIF ! of IF (tracer)
1307
1308	#ifndef MESOSCALE
1309	! update surface pressure
1310	DO ig=1,ngrid
1311	ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
1312	ENDDO
1313
1314	! update pressure levels
1315	DO l=1,nlayer
1316	DO ig=1,ngrid
1317	zplay(ig,l) = aps(l) + bps(l)*ps(ig)
1318	zplev(ig,l) = ap(l) + bp(l)*ps(ig)
1319	ENDDO
1320	ENDDO
1321	zplev(:,nlayer+1) = 0.
1322
1323	! update layers altitude
1324	DO l=2,nlayer
1325	DO ig=1,ngrid
1326	z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
1327	z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
1328	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
1329	ENDDO
1330	ENDDO
1331	#endif
1332
1333	ENDIF ! of IF (callcond)
1334
1335
1336	c-----------------------------------------------------------------------
1337	c 9. Surface and sub-surface soil temperature
1338	c-----------------------------------------------------------------------
1339	c
1340	c
1341	c 9.1 Increment Surface temperature:
1342	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1343
1344	DO ig=1,ngrid
1345	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1346	ENDDO
1347
1348	c Prescribe a cold trap at south pole (except at high obliquity !!)
1349	c Temperature at the surface is set there to be the temperature
1350	c corresponding to equilibrium temperature between phases of CO2
1351
1352
1353	IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
1354	!#ifndef MESOSCALE
1355	! if (caps.and.(obliquit.lt.27.)) then => now done in newcondens
1356	! NB: Updated surface pressure, at grid point 'ngrid', is
1357	! ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1358	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1359	! & (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1360	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095ps(ngrid)))
1361	! endif
1362	!#endif
1363	c -------------------------------------------------------------
1364	c Change of surface albedo in case of ground frost
1365	c everywhere except on the north permanent cap and in regions
1366	c covered by dry ice.
1367	c ALWAYS PLACE these lines after newcondens !!!
1368	c -------------------------------------------------------------
1369	do ig=1,ngrid
1370	if ((co2ice(ig).eq.0).and.
1371	& (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
1372	albedo(ig,1) = albedo_h2o_ice
1373	albedo(ig,2) = albedo_h2o_ice
1374	c write(,) "frost thickness", qsurf(ig,igcm_h2o_ice)
1375	c write(,) "physiq.F frost :"
1376	c & ,lati(ig)180./pi, long(ig)180./pi
1377	endif
1378	enddo ! of do ig=1,ngrid
1379	ENDIF ! of IF (tracer.AND.water.AND.(ngrid.NE.1))
1380
1381	c
1382	c 9.2 Compute soil temperatures and subsurface heat flux:
1383	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1384	IF (callsoil) THEN
1385	c Thermal inertia feedback
1386	IF (tifeedback) THEN
1387	CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
1388	CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
1389	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
1390	ELSE
1391	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1392	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
1393	ENDIF
1394	ENDIF
1395
1396
1397	c-----------------------------------------------------------------------
1398	c 10. Write output files
1399	c ----------------------
1400
1401	c -------------------------------
1402	c Dynamical fields incrementation
1403	c -------------------------------
1404	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1405	! temperature, zonal and meridional wind
1406	DO l=1,nlayer
1407	DO ig=1,ngrid
1408	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1409	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1410	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1411	ENDDO
1412	ENDDO
1413
1414	! tracers
1415	DO iq=1, nq
1416	DO l=1,nlayer
1417	DO ig=1,ngrid
1418	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1419	ENDDO
1420	ENDDO
1421	ENDDO
1422
1423	! Density
1424	DO l=1,nlayer
1425	DO ig=1,ngrid
1426	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1427	ENDDO
1428	ENDDO
1429
1430	! Potential Temperature
1431
1432	DO ig=1,ngrid
1433	DO l=1,nlayer
1434	zh(ig,l) = zt(ig,l)(zplev(ig,1)/zplay(ig,l))*rcp
1435	ENDDO
1436	ENDDO
1437
1438
1439	c Compute surface stress : (NB: z0 is a common in surfdat.h)
1440	c DO ig=1,ngrid
1441	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
1442	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1443	c ENDDO
1444
1445	c Sum of fluxes in solar spectral bands (for output only)
1446	DO ig=1,ngrid
1447	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1448	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1449	ENDDO
1450	c ***** TEST ****************************************************
1451	ztim1 = 999
1452	DO l=1,nlayer
1453	DO ig=1,ngrid
1454	if (pt(ig,l).lt.ztim1) then
1455	ztim1 = pt(ig,l)
1456	igmin = ig
1457	lmin = l
1458	end if
1459	ENDDO
1460	ENDDO
1461	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1462	write(,) 'PHYSIQ: stability WARNING :'
1463	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1464	& 'ig l =', igmin, lmin
1465	end if
1466	c *******************************************************************
1467
1468	c ---------------------
1469	c Outputs to the screen
1470	c ---------------------
1471
1472	IF (lwrite) THEN
1473	PRINT*,'Global diagnostics for the physics'
1474	PRINT,'Variables and their increments x and dx/dt dt'
1475	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1476	WRITE(*,'(2f10.5,2f15.5)')
1477	s tsurf(igout),zdtsurf(igout)*ptimestep,
1478	s zplev(igout,1),pdpsrf(igout)*ptimestep
1479	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1480	WRITE(*,'(i4,6f10.5)') (l,
1481	s pu(igout,l),pdu(igout,l)*ptimestep,
1482	s pv(igout,l),pdv(igout,l)*ptimestep,
1483	s pt(igout,l),pdt(igout,l)*ptimestep,
1484	s l=1,nlayer)
1485	ENDIF ! of IF (lwrite)
1486
1487	c ----------------------------------------------------------
1488	c ----------------------------------------------------------
1489	c INTERPOLATIONS IN THE SURFACE-LAYER
1490	c ----------------------------------------------------------
1491	c ----------------------------------------------------------
1492
1493	n_out=0 ! number of elements in the z_out array.
1494	! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
1495	! to 5
1496	IF (n_out .ne. 0) THEN
1497
1498	IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
1499	IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
1500	IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
1501
1502	z_out(:)=[3.,2.,1.,0.5,0.1]
1503	u_out(:,:)=0.
1504	T_out(:,:)=0.
1505
1506	call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
1507	& g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out,
1508	& T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
1509	! pourquoi ustar recalcule ici? fait dans vdifc.
1510
1511	#ifndef MESOSCALE
1512	IF (ngrid .eq. 1) THEN
1513	dimout=0
1514	ELSE
1515	dimout=2
1516	ENDIF
1517	DO n=1,n_out
1518	write(zstring, '(F8.6)') z_out(n)
1519	call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring),
1520	& 'potential temperature at z_out','K',dimout,T_out(:,n))
1521	call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring),
1522	& 'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n))
1523	ENDDO
1524	call WRITEDIAGFI(ngrid,'u_star',
1525	& 'friction velocity','m/s',dimout,ustar)
1526	call WRITEDIAGFI(ngrid,'teta_star',
1527	& 'friction potential temperature','K',dimout,tstar)
1528	call WRITEDIAGFI(ngrid,'L',
1529	& 'Monin Obukhov length','m',dimout,L_mo)
1530	call WRITEDIAGFI(ngrid,'vvv',
1531	& 'Vertical velocity variance at zout','m',dimout,vvv)
1532	call WRITEDIAGFI(ngrid,'vhf',
1533	& 'Vertical heat flux at zout','m',dimout,vhf)
1534	#else
1535	T_out1(:)=T_out(:,1)
1536	u_out1(:)=u_out(:,1)
1537	#endif
1538
1539	ENDIF
1540
1541	c ----------------------------------------------------------
1542	c ----------------------------------------------------------
1543	c END OF SURFACE LAYER INTERPOLATIONS
1544	c ----------------------------------------------------------
1545	c ----------------------------------------------------------
1546
1547	IF (ngrid.NE.1) THEN
1548
1549	#ifndef MESOSCALE
1550	c -------------------------------------------------------------------
1551	c Writing NetCDF file "RESTARTFI" at the end of the run
1552	c -------------------------------------------------------------------
1553	c Note: 'restartfi' is stored just before dynamics are stored
1554	c in 'restart'. Between now and the writting of 'restart',
1555	c there will have been the itau=itau+1 instruction and
1556	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1557	c thus we store for time=time+dtvr
1558
1559	IF( ((ecritstart.GT.0) .and.
1560	. (MOD(icount*iphysiq,ecritstart).EQ.0))
1561	. .or. lastcall ) THEN
1562
1563	ztime_fin = pday + ptime + ptimestep/(float(iphysiq)*daysec)
1564	. - day_ini - time_phys
1565	print*, pday,ptime,day_ini, time_phys
1566	write(*,'(A,I7,A,F12.5)')
1567	. 'PHYSIQ: Ecriture du fichier restartfi ; icount=',
1568	. icount,' date=',ztime_fin
1569
1570
1571	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
1572	. ptimestep,ztime_fin,
1573	. tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling)
1574
1575	ENDIF
1576	#endif
1577
1578	c -------------------------------------------------------------------
1579	c Calculation of diagnostic variables written in both stats and
1580	c diagfi files
1581	c -------------------------------------------------------------------
1582
1583	if (tracer) then
1584
1585	if(doubleq) then
1586	do ig=1,ngrid
1587	dqdustsurf(ig) =
1588	& zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
1589	dndustsurf(ig) =
1590	& zdqssed(ig,igcm_dust_number)*tauscaling(ig)
1591	ndust(ig,:) =
1592	& pq(ig,:,igcm_dust_number)*tauscaling(ig)
1593	qdust(ig,:) =
1594	& pq(ig,:,igcm_dust_mass)*tauscaling(ig)
1595	enddo
1596	if (scavenging) then
1597	do ig=1,ngrid
1598	dqdustsurf(ig) = dqdustsurf(ig) +
1599	& zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
1600	dndustsurf(ig) = dndustsurf(ig) +
1601	& zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
1602	nccn(ig,:) =
1603	& pq(ig,:,igcm_ccn_number)*tauscaling(ig)
1604	qccn(ig,:) =
1605	& pq(ig,:,igcm_ccn_mass)*tauscaling(ig)
1606	enddo
1607	endif
1608	endif
1609
1610	if (water) then
1611	mtot(:)=0
1612	icetot(:)=0
1613	rave(:)=0
1614	tauTES(:)=0
1615	do ig=1,ngrid
1616	do l=1,nlayer
1617	mtot(ig) = mtot(ig) +
1618	& zq(ig,l,igcm_h2o_vap) *
1619	& (zplev(ig,l) - zplev(ig,l+1)) / g
1620	icetot(ig) = icetot(ig) +
1621	& zq(ig,l,igcm_h2o_ice) *
1622	& (zplev(ig,l) - zplev(ig,l+1)) / g
1623	c Computing abs optical depth at 825 cm-1 in each
1624	c layer to simulate NEW TES retrieval
1625	Qabsice = min(
1626	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
1627	& )
1628	opTES(ig,l)= 0.75 * Qabsice *
1629	& zq(ig,l,igcm_h2o_ice) *
1630	& (zplev(ig,l) - zplev(ig,l+1)) / g
1631	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1632	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1633	enddo
1634	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight
1635	c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1636	enddo
1637	call watersat(ngrid*nlayer,zt,zplay,zqsat)
1638	satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
1639
1640	if (scavenging) then
1641	Nccntot(:)= 0
1642	Mccntot(:)= 0
1643	rave(:)=0
1644	do ig=1,ngrid
1645	do l=1,nlayer
1646	Nccntot(ig) = Nccntot(ig) +
1647	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
1648	& *(zplev(ig,l) - zplev(ig,l+1)) / g
1649	Mccntot(ig) = Mccntot(ig) +
1650	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
1651	& *(zplev(ig,l) - zplev(ig,l+1)) / g
1652	cccc Column integrated effective ice radius
1653	cccc is weighted by total ice surface area (BETTER than total ice mass)
1654	rave(ig) = rave(ig) +
1655	& tauscaling(ig) *
1656	& zq(ig,l,igcm_ccn_number) *
1657	& (zplev(ig,l) - zplev(ig,l+1)) / g *
1658	& rice(ig,l) * rice(ig,l)* (1.+nuice_ref)
1659	enddo
1660	rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
1661	& /max(pi*rave(ig),1.e-30) ! surface weight
1662	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1663	enddo
1664	else ! of if (scavenging)
1665	rave(:)=0
1666	do ig=1,ngrid
1667	do l=1,nlayer
1668	rave(ig) = rave(ig) +
1669	& zq(ig,l,igcm_h2o_ice) *
1670	& (zplev(ig,l) - zplev(ig,l+1)) / g *
1671	& rice(ig,l) * (1.+nuice_ref)
1672	enddo
1673	rave(ig) = max(rave(ig) /
1674	& max(icetot(ig),1.e-30),1.e-30) ! mass weight
1675	enddo
1676	endif ! of if (scavenging)
1677
1678	endif ! of if (water)
1679	endif ! of if (tracer)
1680
1681	#ifndef MESOSCALE
1682	c -----------------------------------------------------------------
1683	c WSTATS: Saving statistics
1684	c -----------------------------------------------------------------
1685	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1686	c which can later be used to make the statistic files of the run:
1687	c "stats") only possible in 3D runs !
1688
1689	IF (callstats) THEN
1690
1691	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1692	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1693	call wstats(ngrid,"co2ice","CO2 ice cover",
1694	& "kg.m-2",2,co2ice)
1695	call wstats(ngrid,"tauref","reference dod at 610 Pa","NU",
1696	& 2,tauref)
1697	call wstats(ngrid,"fluxsurf_lw",
1698	& "Thermal IR radiative flux to surface","W.m-2",2,
1699	& fluxsurf_lw)
1700	call wstats(ngrid,"fluxsurf_sw",
1701	& "Solar radiative flux to surface","W.m-2",2,
1702	& fluxsurf_sw_tot)
1703	call wstats(ngrid,"fluxtop_lw",
1704	& "Thermal IR radiative flux to space","W.m-2",2,
1705	& fluxtop_lw)
1706	call wstats(ngrid,"fluxtop_sw",
1707	& "Solar radiative flux to space","W.m-2",2,
1708	& fluxtop_sw_tot)
1709	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1710	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1711	call wstats(ngrid,"v","Meridional (North-South) wind",
1712	& "m.s-1",3,zv)
1713	call wstats(ngrid,"w","Vertical (down-up) wind",
1714	& "m.s-1",3,pw)
1715	call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
1716	call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
1717	call wstats(ngrid,"q2",
1718	& "Boundary layer eddy kinetic energy",
1719	& "m2.s-2",3,q2)
1720	call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1721	& emis)
1722	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
1723	c & 2,zstress)
1724	c call wstats(ngrid,"sw_htrt","sw heat.rate",
1725	c & "W.m-2",3,zdtsw)
1726	c call wstats(ngrid,"lw_htrt","lw heat.rate",
1727	c & "W.m-2",3,zdtlw)
1728
1729	if (calltherm) then
1730	call wstats(ngrid,"zmax_th","Height of thermals",
1731	& "m",2,zmax_th)
1732	call wstats(ngrid,"hfmax_th","Max thermals heat flux",
1733	& "K.m/s",2,hfmax_th)
1734	call wstats(ngrid,"wstar",
1735	& "Max vertical velocity in thermals",
1736	& "m/s",2,wstar)
1737	endif
1738
1739	if (tracer) then
1740	if (water) then
1741	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
1742	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
1743	call wstats(ngrid,"vmr_h2ovap",
1744	& "H2O vapor volume mixing ratio","mol/mol",
1745	& 3,vmr)
1746	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1747	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
1748	call wstats(ngrid,"vmr_h2oice",
1749	& "H2O ice volume mixing ratio","mol/mol",
1750	& 3,vmr)
1751	vmr=zqsat(1:ngrid,1:nlayer)
1752	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
1753	call wstats(ngrid,"vmr_h2osat",
1754	& "saturation volume mixing ratio","mol/mol",
1755	& 3,vmr)
1756	call wstats(ngrid,"h2o_ice_s",
1757	& "surface h2o_ice","kg/m2",
1758	& 2,qsurf(1,igcm_h2o_ice))
1759	call wstats(ngrid,'albedo',
1760	& 'albedo',
1761	& '',2,albedo(1,1))
1762	call wstats(ngrid,"mtot",
1763	& "total mass of water vapor","kg/m2",
1764	& 2,mtot)
1765	call wstats(ngrid,"icetot",
1766	& "total mass of water ice","kg/m2",
1767	& 2,icetot)
1768	call wstats(ngrid,"reffice",
1769	& "Mean reff","m",
1770	& 2,rave)
1771	call wstats(ngrid,"Nccntot",
1772	& "condensation nuclei","Nbr/m2",
1773	& 2,Nccntot)
1774	call wstats(ngrid,"Mccntot",
1775	& "condensation nuclei mass","kg/m2",
1776	& 2,Mccntot)
1777	call wstats(ngrid,"rice",
1778	& "Ice particle size","m",
1779	& 3,rice)
1780	if (.not.activice) then
1781	call wstats(ngrid,"tauTESap",
1782	& "tau abs 825 cm-1","",
1783	& 2,tauTES)
1784	else
1785	call wstats(ngrid,'tauTES',
1786	& 'tau abs 825 cm-1',
1787	& '',2,taucloudtes)
1788	endif
1789
1790	endif ! of if (water)
1791
1792
1793	if (dustbin.ne.0) then
1794
1795	call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
1796
1797	if (doubleq) then
1798	c call wstats(ngrid,'qsurf','qsurf',
1799	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
1800	c call wstats(ngrid,'Nsurf','N particles',
1801	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
1802	c call wstats(ngrid,'dqsdev','ddevil lift',
1803	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
1804	c call wstats(ngrid,'dqssed','sedimentation',
1805	c & 'kg.m-2.s-1',2,zdqssed(1,1))
1806	c call wstats(ngrid,'dqsdif','diffusion',
1807	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
1808	call wstats(ngrid,'dqsdust',
1809	& 'deposited surface dust mass',
1810	& 'kg.m-2.s-1',2,dqdustsurf)
1811	call wstats(ngrid,'dqndust',
1812	& 'deposited surface dust number',
1813	& 'number.m-2.s-1',2,dndustsurf)
1814	call wstats(ngrid,'reffdust','reffdust',
1815	& 'm',3,rdust*ref_r0)
1816	call wstats(ngrid,'dustq','Dust mass mr',
1817	& 'kg/kg',3,qdust)
1818	call wstats(ngrid,'dustN','Dust number',
1819	& 'part/kg',3,ndust)
1820	else
1821	do iq=1,dustbin
1822	write(str2(1:2),'(i2.2)') iq
1823	call wstats(ngrid,'q'//str2,'mix. ratio',
1824	& 'kg/kg',3,zq(1,1,iq))
1825	call wstats(ngrid,'qsurf'//str2,'qsurf',
1826	& 'kg.m-2',2,qsurf(1,iq))
1827	end do
1828	endif ! (doubleq)
1829
1830	if (scavenging) then
1831	call wstats(ngrid,'ccnq','CCN mass mr',
1832	& 'kg/kg',3,qccn)
1833	call wstats(ngrid,'ccnN','CCN number',
1834	& 'part/kg',3,nccn)
1835	endif ! (scavenging)
1836
1837	endif ! (dustbin.ne.0)
1838
1839
1840	if (thermochem.or.photochem) then
1841	do iq=1,nq
1842	if (noms(iq) .ne. "dust_mass" .and.
1843	$ noms(iq) .ne. "dust_number" .and.
1844	$ noms(iq) .ne. "ccn_mass" .and.
1845	$ noms(iq) .ne. "ccn_number" .and.
1846	$ noms(iq) .ne. "h2o_vap" .and.
1847	$ noms(iq) .ne. "h2o_ice") then
1848	vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
1849	& *mmean(1:ngrid,1:nlayer)/mmol(iq)
1850	rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
1851	& rho(1:ngrid,1:nlayer)n_avog/
1852	& (1000*mmol(iq))
1853	call wstats(ngrid,"vmr_"//trim(noms(iq)),
1854	$ "Volume mixing ratio","mol/mol",3,vmr)
1855	! call wstats(ngrid,"rho_"//trim(noms(iq)),
1856	! $ "Number density","cm-3",3,rhopart)
1857	! call writediagfi(ngrid,"rho_"//trim(noms(iq)),
1858	! $ "Number density","cm-3",3,rhopart)
1859	if ((noms(iq).eq."o") .or. (noms(iq).eq."co2").or.
1860	$ (noms(iq).eq."o3")) then
1861	call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
1862	$ "Volume mixing ratio","mol/mol",3,vmr)
1863	end if
1864	do ig = 1,ngrid
1865	colden(ig,iq) = 0.
1866	end do
1867	do l=1,nlayer
1868	do ig=1,ngrid
1869	colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
1870	$ *(zplev(ig,l)-zplev(ig,l+1))
1871	$ 6.022e22/(mmol(iq)g)
1872	end do
1873	end do
1874	call wstats(ngrid,"c_"//trim(noms(iq)),
1875	$ "column","mol cm-2",2,colden(1,iq))
1876	call writediagfi(ngrid,"c_"//trim(noms(iq)),
1877	$ "column","mol cm-2",2,colden(1,iq))
1878	end if ! of if (noms(iq) .ne. "dust_mass" ...)
1879	end do ! of do iq=1,nq
1880	end if ! of if (thermochem.or.photochem)
1881
1882	end if ! of if (tracer)
1883
1884	IF(lastcall) THEN
1885	write (,) "Writing stats..."
1886	call mkstats(ierr)
1887	ENDIF
1888
1889	ENDIF !if callstats
1890
1891	c (Store EOF for Mars Climate database software)
1892	IF (calleofdump) THEN
1893	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1894	ENDIF
1895	#endif
1896	!endif of ifndef MESOSCALE
1897
1898	#ifdef MESOSCALE
1899
1900	!! see comm_wrf.
1901	!! not needed when an array is already in a shared module.
1902	!! --> example : hfmax_th, zmax_th
1903
1904	!state real HR_SW ikj misc 1 - h "HR_SW" "HEATING RATE SW" "K/s"
1905	comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
1906	!state real HR_LW ikj misc 1 - h "HR_LW" "HEATING RATE LW" "K/s"
1907	comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
1908	!state real SWDOWNZ ij misc 1 - h "SWDOWNZ" "DOWNWARD SW FLUX AT SURFACE" "W m-2"
1909	comm_SWDOWNZ(1:ngrid) = fluxsurf_sw_tot(1:ngrid)
1910	!state real TAU_DUST ij misc 1 - h "TAU_DUST" "REFERENCE VISIBLE DUST OPACITY" ""
1911	comm_TAU_DUST(1:ngrid) = tauref(1:ngrid)
1912	!state real RDUST ikj misc 1 - h "RDUST" "DUST RADIUS" "m"
1913	comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer)
1914	!state real QSURFDUST ij misc 1 - h "QSURFDUST" "DUST MASS AT SURFACE" "kg m-2"
1915	IF (igcm_dust_mass .ne. 0) THEN
1916	comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
1917	ELSE
1918	comm_QSURFDUST(1:ngrid) = 0.
1919	ENDIF
1920	!state real MTOT ij misc 1 - h "MTOT" "TOTAL MASS WATER VAPOR in pmic" "pmic"
1921	comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
1922	!state real ICETOT ij misc 1 - h "ICETOT" "TOTAL MASS WATER ICE" "kg m-2"
1923	comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
1924	!state real VMR_ICE ikj misc 1 - h "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm"
1925	IF (igcm_h2o_ice .ne. 0) THEN
1926	comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6
1927	. * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1928	. * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice)
1929	ELSE
1930	comm_VMR_ICE(1:ngrid,1:nlayer) = 0.
1931	ENDIF
1932	!state real TAU_ICE ij misc 1 - h "TAU_ICE" "CLOUD OD at 825 cm-1 TES" ""
1933	comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
1934	!state real RICE ikj misc 1 - h "RICE" "ICE RADIUS" "m"
1935	comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)
1936
1937
1938	!! calculate sensible heat flux in W/m2 for outputs
1939	!! -- the one computed in vdifc is not the real one
1940	!! -- vdifc must have been called
1941	if (.not.callrichsl) then
1942	sensibFlux(1:ngrid) = zflubid(1:ngrid)
1943	. - capcal(1:ngrid)*zdtsdif(1:ngrid)
1944	else
1945	sensibFlux(1:ngrid) =
1946	& (pplay(1:ngrid,1)/(rpt(1:ngrid,1)))cpp
1947	& sqrt(pu(1:ngrid,1)pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1)
1948	& +(log(1.+0.7wstar(1:ngrid) + 2.3wstar(1:ngrid)2))2)
1949	& zcdh(1:ngrid)(tsurf(1:ngrid)-zh(1:ngrid,1))
1950	endif
1951
1952	#else
1953	#ifndef MESOINI
1954
1955	c ==========================================================
1956	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
1957	c any variable for diagnostic (output with period
1958	c "ecritphy", set in "run.def")
1959	c ==========================================================
1960	c WRITEDIAGFI can ALSO be called from any other subroutines
1961	c for any variables !!
1962	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1963	c & emis)
1964	c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
1965	c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
1966	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1967	& tsurf)
1968	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1969	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
1970	& ,"kg.m-2",2,co2ice)
1971
1972	call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
1973	& "K",2,zt(1,7))
1974	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1975	& fluxsurf_lw)
1976	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1977	& fluxsurf_sw_tot)
1978	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1979	& fluxtop_lw)
1980	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1981	& fluxtop_sw_tot)
1982	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1983	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1984	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1985	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1986	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1987	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1988	c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1989	call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay)
1990	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
1991	c & zstress)
1992	c call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
1993	c & 'w.m-2',3,zdtsw)
1994	c call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
1995	c & 'w.m-2',3,zdtlw)
1996	if (.not.activice) then
1997	CALL WRITEDIAGFI(ngrid,'tauTESap',
1998	& 'tau abs 825 cm-1',
1999	& '',2,tauTES)
2000	else
2001	CALL WRITEDIAGFI(ngrid,'tauTES',
2002	& 'tau abs 825 cm-1',
2003	& '',2,taucloudtes)
2004	endif
2005
2006	#else
2007	!!! this is to ensure correct initialisation of mesoscale model
2008	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
2009	& tsurf)
2010	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
2011	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
2012	& co2ice)
2013	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
2014	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
2015	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
2016	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
2017	& emis)
2018	call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
2019	& "K",3,tsoil)
2020	call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
2021	& "K",3,inertiedat)
2022	#endif
2023
2024
2025	c ----------------------------------------------------------
2026	c Outputs of the CO2 cycle
2027	c ----------------------------------------------------------
2028
2029	if (tracer.and.(igcm_co2.ne.0)) then
2030	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
2031	! & "kg/kg",2,zq(1,1,igcm_co2))
2032	call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
2033	& "kg/kg",3,zq(1,1,igcm_co2))
2034
2035	! Compute co2 column
2036	co2col(:)=0
2037	do l=1,nlayer
2038	do ig=1,ngrid
2039	co2col(ig)=co2col(ig)+
2040	& zq(ig,l,igcm_co2)*(zplev(ig,l)-zplev(ig,l+1))/g
2041	enddo
2042	enddo
2043	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
2044	& co2col)
2045	endif ! of if (tracer.and.(igcm_co2.ne.0))
2046
2047	c ----------------------------------------------------------
2048	c Outputs of the water cycle
2049	c ----------------------------------------------------------
2050	if (tracer) then
2051	if (water) then
2052
2053	#ifdef MESOINI
2054	!!!! waterice = q01, voir readmeteo.F90
2055	call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
2056	& 'kg/kg',3,
2057	& zq(1:ngrid,1:nlayer,igcm_h2o_ice))
2058	!!!! watervapor = q02, voir readmeteo.F90
2059	call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
2060	& 'kg/kg',3,
2061	& zq(1:ngrid,1:nlayer,igcm_h2o_vap))
2062	!!!! surface waterice qsurf02 (voir readmeteo)
2063	call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
2064	& 'kg.m-2',2,
2065	& qsurf(1:ngrid,igcm_h2o_ice))
2066	#endif
2067
2068	CALL WRITEDIAGFI(ngrid,'mtot',
2069	& 'total mass of water vapor',
2070	& 'kg/m2',2,mtot)
2071	CALL WRITEDIAGFI(ngrid,'icetot',
2072	& 'total mass of water ice',
2073	& 'kg/m2',2,icetot)
2074	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2075	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
2076	call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
2077	& 'mol/mol',3,vmr)
2078	vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
2079	& *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2080	call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
2081	& 'mol/mol',3,vmr)
2082	CALL WRITEDIAGFI(ngrid,'reffice',
2083	& 'Mean reff',
2084	& 'm',2,rave)
2085	if (scavenging) then
2086	CALL WRITEDIAGFI(ngrid,"Nccntot",
2087	& "condensation nuclei","Nbr/m2",
2088	& 2,Nccntot)
2089	CALL WRITEDIAGFI(ngrid,"Mccntot",
2090	& "mass condensation nuclei","kg/m2",
2091	& 2,Mccntot)
2092	endif
2093	call WRITEDIAGFI(ngrid,'rice','Ice particle size',
2094	& 'm',3,rice)
2095	call WRITEDIAGFI(ngrid,'h2o_ice_s',
2096	& 'surface h2o_ice',
2097	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
2098	CALL WRITEDIAGFI(ngrid,'albedo',
2099	& 'albedo',
2100	& '',2,albedo(1,1))
2101	if (tifeedback) then
2102	call WRITEDIAGSOIL(ngrid,"soiltemp",
2103	& "Soil temperature","K",
2104	& 3,tsoil)
2105	call WRITEDIAGSOIL(ngrid,'soilti',
2106	& 'Soil Thermal Inertia',
2107	& 'J.s-1/2.m-2.K-1',3,inertiesoil)
2108	endif
2109	endif !(water)
2110
2111
2112	if (water.and..not.photochem) then
2113	iq=nq
2114	c write(str2(1:2),'(i2.2)') iq
2115	c call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
2116	c & 'kg.m-2',2,zdqscloud(1,iq))
2117	c call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
2118	c & 'kg/kg',3,zdqchim(1,1,iq))
2119	c call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
2120	c & 'kg/kg',3,zdqdif(1,1,iq))
2121	c call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
2122	c & 'kg/kg',3,zdqadj(1,1,iq))
2123	c call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
2124	c & 'kg/kg',3,zdqc(1,1,iq))
2125	endif !(water.and..not.photochem)
2126	endif
2127
2128	c ----------------------------------------------------------
2129	c Outputs of the dust cycle
2130	c ----------------------------------------------------------
2131
2132	call WRITEDIAGFI(ngrid,'tauref',
2133	& 'Dust ref opt depth','NU',2,tauref)
2134
2135	if (tracer.and.(dustbin.ne.0)) then
2136	call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
2137	if (doubleq) then
2138	c call WRITEDIAGFI(ngrid,'qsurf','qsurf',
2139	c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))
2140	c call WRITEDIAGFI(ngrid,'Nsurf','N particles',
2141	c & 'N.m-2',2,qsurf(1,igcm_dust_number))
2142	c call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift',
2143	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
2144	c call WRITEDIAGFI(ngrid,'dqssed','sedimentation',
2145	c & 'kg.m-2.s-1',2,zdqssed(1,1))
2146	c call WRITEDIAGFI(ngrid,'dqsdif','diffusion',
2147	c & 'kg.m-2.s-1',2,zdqsdif(1,1))
2148	c call WRITEDIAGFI(ngrid,'sedice','sedimented ice',
2149	c & 'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice))
2150	c call WRITEDIAGFI(ngrid,'subice','sublimated ice',
2151	c & 'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice))
2152	call WRITEDIAGFI(ngrid,'dqsdust',
2153	& 'deposited surface dust mass',
2154	& 'kg.m-2.s-1',2,dqdustsurf)
2155	call WRITEDIAGFI(ngrid,'dqndust',
2156	& 'deposited surface dust number',
2157	& 'number.m-2.s-1',2,dndustsurf)
2158	call WRITEDIAGFI(ngrid,'reffdust','reffdust',
2159	& 'm',3,rdust*ref_r0)
2160	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2161	& 'kg/kg',3,qdust)
2162	call WRITEDIAGFI(ngrid,'dustN','Dust number',
2163	& 'part/kg',3,ndust)
2164	c call WRITEDIAGFI(ngrid,"tauscaling",
2165	c & "dust conversion factor"," ",2,tauscaling)
2166	#ifdef MESOINI
2167	! !!! to initialize mesoscale we need scaled variables
2168	! !!! because this must correspond to starting point for tracers
2169	! call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2170	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
2171	! call WRITEDIAGFI(ngrid,'dustN','Dust number',
2172	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
2173	! call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
2174	! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
2175	! call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
2176	! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
2177	call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2178	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
2179	call WRITEDIAGFI(ngrid,'dustN','Dust number',
2180	& 'part/kg',3,pq(1,1,igcm_dust_number))
2181	call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
2182	& 'kg/kg',3,pq(1,1,igcm_ccn_mass))
2183	call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
2184	& 'part/kg',3,pq(1,1,igcm_ccn_number))
2185	#endif
2186	else
2187	do iq=1,dustbin
2188	write(str2(1:2),'(i2.2)') iq
2189	call WRITEDIAGFI(ngrid,'q'//str2,'mix. ratio',
2190	& 'kg/kg',3,zq(1,1,iq))
2191	call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
2192	& 'kg.m-2',2,qsurf(1,iq))
2193	end do
2194	endif ! (doubleq)
2195
2196	if (scavenging) then
2197	call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
2198	& 'kg/kg',3,qccn)
2199	call WRITEDIAGFI(ngrid,'ccnN','CCN number',
2200	& 'part/kg',3,nccn)
2201	endif ! (scavenging)
2202
2203	c if (submicron) then
2204	c call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',
2205	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
2206	c endif ! (submicron)
2207	end if ! (tracer.and.(dustbin.ne.0))
2208
2209
2210	c ----------------------------------------------------------
2211	c Thermospheric outputs
2212	c ----------------------------------------------------------
2213
2214	if(callthermos) then
2215
2216	call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
2217	$ 3,zdtnlte)
2218	call WRITEDIAGFI(ngrid,"quv","UV heating","K/s",
2219	$ 3,zdteuv)
2220	call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s",
2221	$ 3,zdtconduc)
2222	call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
2223	$ 3,zdtnirco2)
2224
2225	endif !(callthermos)
2226
2227	c ----------------------------------------------------------
2228	c ----------------------------------------------------------
2229	c PBL OUTPUS
2230	c ----------------------------------------------------------
2231	c ----------------------------------------------------------
2232
2233	c ----------------------------------------------------------
2234	c Outputs of thermals
2235	c ----------------------------------------------------------
2236	if (calltherm) then
2237
2238	! call WRITEDIAGFI(ngrid,'dtke',
2239	! & 'tendance tke thermiques','m2/s2',
2240	! & 3,dtke_th)
2241	! call WRITEDIAGFI(ngrid,'d_u_ajs',
2242	! & 'tendance u thermiques','m/s',
2243	! & 3,pdu_th*ptimestep)
2244	! call WRITEDIAGFI(ngrid,'d_v_ajs',
2245	! & 'tendance v thermiques','m/s',
2246	! & 3,pdv_th*ptimestep)
2247	! if (tracer) then
2248	! if (nq .eq. 2) then
2249	! call WRITEDIAGFI(ngrid,'deltaq_th',
2250	! & 'delta q thermiques','kg/kg',
2251	! & 3,ptimestep*pdq_th(:,:,2))
2252	! endif
2253	! endif
2254
2255	call WRITEDIAGFI(ngrid,'zmax_th',
2256	& 'hauteur du thermique','m',
2257	& 2,zmax_th)
2258	call WRITEDIAGFI(ngrid,'hfmax_th',
2259	& 'maximum TH heat flux','K.m/s',
2260	& 2,hfmax_th)
2261	call WRITEDIAGFI(ngrid,'wstar',
2262	& 'maximum TH vertical velocity','m/s',
2263	& 2,wstar)
2264
2265	endif
2266
2267	c ----------------------------------------------------------
2268	c ----------------------------------------------------------
2269	c END OF PBL OUTPUS
2270	c ----------------------------------------------------------
2271	c ----------------------------------------------------------
2272
2273
2274	c ----------------------------------------------------------
2275	c Output in netcdf file "diagsoil.nc" for subterranean
2276	c variables (output every "ecritphy", as for writediagfi)
2277	c ----------------------------------------------------------
2278
2279	! Write soil temperature
2280	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
2281	! & 3,tsoil)
2282	! Write surface temperature
2283	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
2284	! & 2,tsurf)
2285
2286	c ==========================================================
2287	c END OF WRITEDIAGFI
2288	c ==========================================================
2289	#endif
2290	! of ifdef MESOSCALE
2291
2292	ELSE ! if(ngrid.eq.1)
2293
2294	#ifndef MESOSCALE
2295	write(*,'("Ls =",f11.6," tauref(",f4.0," Pa) =",f9.6)')
2296	& zls*180./pi,odpref,tauref
2297	c ----------------------------------------------------------------------
2298	c Output in grads file "g1d" (ONLY when using testphys1d)
2299	c (output at every X physical timestep)
2300	c ----------------------------------------------------------------------
2301	c
2302	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
2303	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
2304	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
2305
2306	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
2307	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
2308	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
2309	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
2310
2311	! THERMALS STUFF 1D
2312	if(calltherm) then
2313
2314	call WRITEDIAGFI(ngrid,'lmax_th',
2315	& 'hauteur du thermique','point',
2316	& 0,lmax_th_out)
2317	call WRITEDIAGFI(ngrid,'zmax_th',
2318	& 'hauteur du thermique','m',
2319	& 0,zmax_th)
2320	call WRITEDIAGFI(ngrid,'hfmax_th',
2321	& 'maximum TH heat flux','K.m/s',
2322	& 0,hfmax_th)
2323	call WRITEDIAGFI(ngrid,'wstar',
2324	& 'maximum TH vertical velocity','m/s',
2325	& 0,wstar)
2326
2327	co2col(:)=0.
2328	if (tracer) then
2329	do l=1,nlayer
2330	do ig=1,ngrid
2331	co2col(ig)=co2col(ig)+
2332	& zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
2333	enddo
2334	enddo
2335
2336	end if
2337	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
2338	& ,'kg/m-2',0,co2col)
2339	endif ! of if (calltherm)
2340
2341	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
2342	& ,'m/s',1,pw)
2343	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
2344	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
2345	& tsurf)
2346	call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
2347	call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
2348
2349	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
2350	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
2351	call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
2352	! call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate", &
2353	! & "K.s-1",1,dtrad/zpopsk)
2354	! call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
2355	! & 'w.m-2',1,zdtsw/zpopsk)
2356	! call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
2357	! & 'w.m-2',1,zdtlw/zpopsk)
2358	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
2359	& ,"kg.m-2",0,co2ice)
2360
2361	! or output in diagfi.nc (for testphys1d)
2362	call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
2363	call WRITEDIAGFI(ngrid,'temp','Temperature',
2364	& 'K',1,zt)
2365
2366	if(tracer) then
2367	c CALL writeg1d(ngrid,1,tau,'tau','SI')
2368	do iq=1,nq
2369	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
2370	call WRITEDIAGFI(ngrid,trim(noms(iq)),
2371	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
2372	end do
2373	if (doubleq) then
2374	call WRITEDIAGFI(ngrid,'rdust','rdust',
2375	& 'm',1,rdust)
2376	endif
2377	if (water.AND.tifeedback) then
2378	call WRITEDIAGFI(ngrid,"soiltemp",
2379	& "Soil temperature","K",
2380	& 1,tsoil)
2381	call WRITEDIAGFI(ngrid,'soilti',
2382	& 'Soil Thermal Inertia',
2383	& 'J.s-1/2.m-2.K-1',1,inertiesoil)
2384	endif
2385	end if
2386
2387	cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
2388	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2389	IF (water) THEN
2390
2391	if (.not.activice) then
2392
2393	tauTES=0
2394	do l=1,nlayer
2395	Qabsice = min(
2396	& max(0.4e6rice(1,l)(1.+nuice_ref)-0.05 ,0.),1.2
2397	& )
2398	opTES(1,l)= 0.75 * Qabsice *
2399	& zq(1,l,igcm_h2o_ice) *
2400	& (zplev(1,l) - zplev(1,l+1)) / g
2401	& / (rho_ice * rice(1,l) * (1.+nuice_ref))
2402	tauTES=tauTES+ opTES(1,l)
2403	enddo
2404	CALL WRITEDIAGFI(ngrid,'tauTESap',
2405	& 'tau abs 825 cm-1',
2406	& '',0,tauTES)
2407	else
2408
2409	CALL WRITEDIAGFI(ngrid,'tauTES',
2410	& 'tau abs 825 cm-1',
2411	& '',0,taucloudtes)
2412	endif
2413
2414	mtot = 0
2415	icetot = 0
2416	h2otot = qsurf(1,igcm_h2o_ice)
2417
2418	do l=1,nlayer
2419	mtot = mtot + zq(1,l,igcm_h2o_vap)
2420	& * (zplev(1,l) - zplev(1,l+1)) / g
2421	icetot = icetot + zq(1,l,igcm_h2o_ice)
2422	& * (zplev(1,l) - zplev(1,l+1)) / g
2423	end do
2424	h2otot = h2otot+mtot+icetot
2425
2426	CALL WRITEDIAGFI(ngrid,'h2otot',
2427	& 'h2otot',
2428	& 'kg/m2',0,h2otot)
2429	CALL WRITEDIAGFI(ngrid,'mtot',
2430	& 'mtot',
2431	& 'kg/m2',0,mtot)
2432	CALL WRITEDIAGFI(ngrid,'icetot',
2433	& 'icetot',
2434	& 'kg/m2',0,icetot)
2435
2436	if (scavenging) then
2437
2438	rave = 0
2439	do l=1,nlayer
2440	cccc Column integrated effective ice radius
2441	cccc is weighted by total ice surface area (BETTER)
2442	rave = rave + tauscaling(1) *
2443	& zq(1,l,igcm_ccn_number) *
2444	& (zplev(1,l) - zplev(1,l+1)) / g *
2445	& rice(1,l) * rice(1,l)* (1.+nuice_ref)
2446	enddo
2447	rave=icetot0.75/max(ravepi*rho_ice,1.e-30) ! surface weight
2448
2449	Nccntot= 0
2450	call watersat(ngrid*nlayer,zt,zplay,zqsat)
2451	do l=1,nlayer
2452	Nccntot = Nccntot +
2453	& zq(1,l,igcm_ccn_number)*tauscaling(1)
2454	& *(zplev(1,l) - zplev(1,l+1)) / g
2455	satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
2456	satu(1,l) = (max(satu(1,l),float(1))-1)
2457	! & * zq(1,l,igcm_h2o_vap) *
2458	! & (zplev(1,l) - zplev(1,l+1)) / g
2459	enddo
2460	call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
2461	& satu)
2462	CALL WRITEDIAGFI(ngrid,'Nccntot',
2463	& 'Nccntot',
2464	& 'nbr/m2',0,Nccntot)
2465
2466	call WRITEDIAGFI(ngrid,'zdqsed_dustq'
2467	& ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
2468	call WRITEDIAGFI(ngrid,'zdqsed_dustN'
2469	&,'sedimentation N','Nbr.m-2.s-1',1,
2470	& zdqsed(1,:,igcm_dust_number))
2471
2472	else ! of if (scavenging)
2473
2474	cccc Column integrated effective ice radius
2475	cccc is weighted by total ice mass (LESS GOOD)
2476	rave = 0
2477	do l=1,nlayer
2478	rave = rave + zq(1,l,igcm_h2o_ice)
2479	& * (zplev(1,l) - zplev(1,l+1)) / g
2480	& * rice(1,l) * (1.+nuice_ref)
2481	enddo
2482	rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight
2483	endif ! of if (scavenging)
2484
2485	CALL WRITEDIAGFI(ngrid,'reffice',
2486	& 'reffice',
2487	& 'm',0,rave)
2488
2489	do iq=1,nq
2490	call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
2491	& trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
2492	end do
2493
2494	call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
2495	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
2496	call WRITEDIAGFI(ngrid,'zdqcloud_vap','cloud vap',
2497	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
2498	call WRITEDIAGFI(ngrid,'zdqcloud','cloud ice',
2499	& 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
2500	& +zdqcloud(1,:,igcm_h2o_vap))
2501
2502	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
2503	& rice)
2504	ENDIF ! of IF (water)
2505
2506	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2507	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2508
2509
2510	zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
2511
2512	do l=2,nlayer-1
2513	tmean=zt(1,l)
2514	if(zt(1,l).ne.zt(1,l-1))
2515	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
2516	zlocal(l)= zlocal(l-1)
2517	& -log(zplay(1,l)/zplay(1,l-1))rnew(1,l)tmean/g
2518	enddo
2519	zlocal(nlayer)= zlocal(nlayer-1)-
2520	& log(zplay(1,nlayer)/zplay(1,nlayer-1))*
2521	& rnew(1,nlayer)*tmean/g
2522	#endif
2523
2524	END IF ! if(ngrid.ne.1)
2525
2526	icount=icount+1
2527	RETURN
2528	END

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

Download in other formats:

Original Format