Context Navigation

physiq.F @ 636

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

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

Download in other formats:

Original Format