Context Navigation

physiq.F @ 461

Last change on this file since 461 was 459, checked in by emillour, 14 years ago

Mars GCM: Include changes and updates for photochemistry by FL:

aeronomars/calchim.F : change in units of surface density.
aeronomars/surfacearea.F : new routine to compute ice and dust surface area

(m2/m3) available for heterogeneous reactions.

phymars/initracer.F : bug correction: initialize igcm_ch4 and change loop

bounds (when guessing tracer names/properties with
old input files).

phymars/watercloud.F : cleanup.
phymars/physiq.F : add call to surfacearea; photochemistry is now called

after sedimentation to take into acount updated rdust
and rice.

File size: 74.2 KB

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

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

Download in other formats:

Original Format