Context Navigation

physiq.F @ 414

Last change on this file since 414 was 414, checked in by aslmd, 14 years ago

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F

File size: 74.3 KB

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

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

Download in other formats:

Original Format