Context Navigation

physiq.F @ 496

Last change on this file since 496 was 496, checked in by acolaitis, 13 years ago

M 489 libf/phymars/thermcell_main_mars.F90
^{--------------------> New parameters for thermals, following latest version of the relevant article}

M 489 libf/phymars/physiq.F
^{--------------------> Minor modifications}

D 489 libf/phymars/surflayer_interpol.F
A 0 libf/phymars/pbl_parameters.F
^{--------------------> Replaced surflayer_interpol.F by a cleaner and richer version called pbl_parameters.F}

This routine is the base of what will later be implemented in the MCD and as a tool to
compute surface properties from output fields (so that it will ultimately disappear from libf
and might become an utility)

M 489 libf/phymars/vdif_cd.F
^{--------------------> Minor modification to the Richardson computation}

M 489 libf/phymars/vdifc.F
^{--------------------> Minor modification to the aerodynamic conductances computation, replacing wmax by hfmax, which}

is a better proxy for convective activity. Might be replaced by w_star later.

File size: 73.5 KB

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

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

Download in other formats:

Original Format