Context Navigation

physiq.F @ 532

Last change on this file since 532 was 529, checked in by acolaitis, 13 years ago
Physiq: changed called to pbl_parameter and vdifc. Vdifc: added a flag to deactivate mass-variation scheme for LES computations for two reason: it is costy for nothing and does not work in polar LES.
File size: 77.7 KB

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

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

Download in other formats:

Original Format