Context Navigation

physiq.F @ 1266

Last change on this file since 1266 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

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

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

Download in other formats:

Original Format