Context Navigation

physiq.F @ 563

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

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

Download in other formats:

Original Format