Context Navigation

physiq.F @ 520

Last change on this file since 520 was 520, checked in by tnavarro, 13 years ago

10/02/12 == TN

Major update on watercycle: a smaller integration timestep is now used

in watercloud.F, sedimentation of clouds is done in watercloud instead of
callsedim.F

Temperature-dependant contact parameter in nuclea.F
No dust lifting if CO2 ice in vdif.c
Ice integrated column opacity is written in diagfi from physiq.F, instead

of aeropacity.F. Mandatory if iradia is not 1.

New definition of permanent ice in surfini.F and possibility to have an ice

cap in it in 1d.

Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging

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

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

Download in other formats:

Original Format