Context Navigation

meso_physiq.F @ 227

Last change on this file since 227 was 227, checked in by aslmd, 13 years ago
MESOSCALE and LMDZ.MARS: corrected the bug with readtesassim, now the GCM routine can be used. changed name for slope routines to drop the meso_ prefix
File size: 59.3 KB

Line
1	SUBROUTINE meso_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 version: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
67	c
68	c arguments:
69	c ----------
70	c
71	c input:
72	c ------
73	c ecri period (in dynamical timestep) to write output
74	c ngrid Size of the horizontal grid.
75	c All internal loops are performed on that grid.
76	c nlayer Number of vertical layers.
77	c nq Number of advected fields
78	c firstcall True at the first call
79	c lastcall True at the last call
80	c pday Number of days counted from the North. Spring
81	c equinoxe.
82	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
83	c ptimestep timestep (s)
84	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
85	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
86	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
87	c pu(ngrid,nlayer) u component of the wind (ms-1)
88	c pv(ngrid,nlayer) v component of the wind (ms-1)
89	c pt(ngrid,nlayer) Temperature (K)
90	c pq(ngrid,nlayer,nq) Advected fields
91	c pudyn(ngrid,nlayer) \
92	c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
93	c ptdyn(ngrid,nlayer) / corresponding variables
94	c pqdyn(ngrid,nlayer,nq) /
95	c pw(ngrid,?) vertical velocity
96	c
97	c output:
98	c -------
99	c
100	c pdu(ngrid,nlayermx) \
101	c pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
102	c pdt(ngrid,nlayermx) / variables due to physical processes.
103	c pdq(ngrid,nlayermx,nqmx) /
104	c pdpsrf(ngrid) /
105	c tracerdyn call tracer in dynamical part of GCM ?
106
107	c
108	c=======================================================================
109	c
110	c 0. Declarations :
111	c ------------------
112
113	#include "dimensions.h"
114	#include "dimphys.h"
115	#include "comgeomfi.h"
116	#include "surfdat.h"
117	#include "comsoil.h"
118	#include "comdiurn.h"
119	#include "callkeys.h"
120	#include "comcstfi.h"
121	#include "planete.h"
122	#include "comsaison.h"
123	#include "control.h"
124	#include "dimradmars.h"
125	#include "comg1d.h"
126	#include "tracer.h"
127	#include "nlteparams.h"
128
129	#include "chimiedata.h"
130	#include "watercap.h"
131	#include "param.h"
132	#include "param_v3.h"
133	#include "conc.h"
134
135	#include "netcdf.inc"
136
137	#ifdef MESOSCALE
138	#include "slope.h"
139	#include "wrf_output_2d.h"
140	#include "wrf_output_3d.h"
141	#include "advtrac.h" !!! this is necessary for tracers (in dyn3d)
142	#include "meso_inc/meso_inc_var.F"
143	#endif
144
145	c Arguments :
146	c -----------
147
148	c inputs:
149	c -------
150	INTEGER ngrid,nlayer,nq
151	REAL ptimestep
152	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
153	REAL pphi(ngridmx,nlayer)
154	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
155	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
156	REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
157	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
158	LOGICAL firstcall,lastcall
159
160	REAL pday
161	REAL ptime
162	logical tracerdyn
163
164	c outputs:
165	c --------
166	c physical tendencies
167	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
168	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
169	REAL pdpsrf(ngridmx) ! surface pressure tendency
170
171
172	c Local saved variables:
173	c ----------------------
174	c aerosol (dust or ice) extinction optical depth at reference wavelength
175	c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
176	c aerosol optical properties :
177	REAL aerosol(ngridmx,nlayermx,naerkind)
178
179	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
180	INTEGER icount ! counter of calls to physiq during the run.
181	REAL tsurf(ngridmx) ! Surface temperature (K)
182	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
183	REAL co2ice(ngridmx) ! co2 ice surface layer (kg.m-2)
184	REAL albedo(ngridmx,2) ! Surface albedo in each solar band
185	REAL emis(ngridmx) ! Thermal IR surface emissivity
186	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
187	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
188	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
189	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
190	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
191	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
192	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
193	INTEGER ig_vl1 ! Grid Point near VL1 (for diagnostic)
194
195	c Variables used by the water ice microphysical scheme:
196	REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m)
197	REAL nuice(ngridmx,nlayermx) ! Estimated effective variance
198	! of the size distribution
199	c Albedo of deposited surface ice
200	!!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
201	REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB
202
203	SAVE day_ini, icount
204	SAVE aerosol, tsurf,tsoil
205	SAVE co2ice,albedo,emis, q2
206	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
207	SAVE ig_vl1
208
209	REAL stephan
210	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
211	SAVE stephan
212
213	c Local variables :
214	c -----------------
215
216	REAL CBRT
217	EXTERNAL CBRT
218
219	CHARACTER*80 fichier
220	INTEGER l,ig,ierr,igout,iq,i, tapphys
221
222	REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)
223	REAL fluxsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)
224	REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)
225	REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)
226	REAL tauref(ngridmx) ! Reference column optical depth at 700 Pa
227	! (used if active=F)
228	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
229	REAL zls ! solar longitude (rad)
230	REAL zday ! date (time since Ls=0, in martian days)
231	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
232	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
233	REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
234
235	c Tendancies due to various processes:
236	REAL dqsurf(ngridmx,nqmx)
237	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
238	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
239	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear area
240	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear area
241	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
242	REAL zdtnlte(ngridmx,nlayermx) ! (K/s)
243	REAL zdtsurf(ngridmx) ! (K/s)
244	REAL zdtcloud(ngridmx,nlayermx)
245	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
246	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
247	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
248	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
249	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
250	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
251	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
252	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
253
254	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
255	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
256	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
257	REAL zdqadj(ngridmx,nlayermx,nqmx)
258	REAL zdqc(ngridmx,nlayermx,nqmx)
259	REAL zdqcloud(ngridmx,nlayermx,nqmx)
260	REAL zdqscloud(ngridmx,nqmx)
261	REAL zdqchim(ngridmx,nlayermx,nqmx)
262	REAL zdqschim(ngridmx,nqmx)
263
264	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
265	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
266	REAL zdumolvis(ngridmx,nlayermx)
267	REAL zdvmolvis(ngridmx,nlayermx)
268	real zdqmoldiff(ngridmx,nlayermx,nqmx)
269
270	c Local variable for local intermediate calcul:
271	REAL zflubid(ngridmx)
272	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
273	REAL zdum1(ngridmx,nlayermx)
274	REAL zdum2(ngridmx,nlayermx)
275	REAL ztim1,ztim2,ztim3, z1,z2
276	REAL ztime_fin
277	REAL zdh(ngridmx,nlayermx)
278	INTEGER length
279	PARAMETER (length=100)
280
281	c local variables only used for diagnostic (output in file "diagfi" or "stats")
282	c -----------------------------------------------------------------------------
283	REAL ps(ngridmx), zt(ngridmx,nlayermx)
284	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
285	REAL zq(ngridmx,nlayermx,nqmx)
286	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
287	character*2 str2
288	character*5 str5
289	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
290	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
291	! (particules kg-1)
292	SAVE ccn !! in case iradia != 1
293	real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
294	real qtot1,qtot2 ! total aerosol mass
295	integer igmin, lmin
296	logical tdiag
297
298	real co2col(ngridmx) ! CO2 column
299	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
300	REAL zstress(ngrid), cd
301	real hco2(nqmx),tmean, zlocal(nlayermx)
302	real rho(ngridmx,nlayermx) ! density
303	real vmr(ngridmx,nlayermx) ! volume mixing ratio
304	REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2)
305	REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
306	REAL rave(ngridmx) ! Mean water ice effective radius (m)
307	REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
308	REAL tauTES(ngridmx) ! column optical depth at 825 cm-1
309	REAL Qabsice ! Water ice absorption coefficient
310
311
312	REAL time_phys
313
314	c Variables from thermal
315
316	REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
317	REAL wmax_th(ngridmx)
318	REAL ,SAVE :: hfmax_th(ngridmx)
319	REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
320	REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
321	INTEGER lmax_th(ngridmx)
322	REAL dtke_th(ngridmx,nlayermx+1)
323	REAL dummycol(ngridmx)
324
325	c=======================================================================
326
327	c 1. Initialisation:
328	c -----------------
329
330	c 1.1 Initialisation only at first call
331	c ---------------------------------------
332	IF (firstcall) THEN
333
334	c variables set to 0
335	c ~~~~~~~~~~~~~~~~~~
336	call zerophys(ngridnlayernaerkind,aerosol)
337	call zerophys(ngrid*nlayer,dtrad)
338	call zerophys(ngrid,fluxrad)
339
340	c read startfi
341	c ~~~~~~~~~~~~
342	#ifndef MESOSCALE
343	! Read netcdf initial physical parameters.
344	CALL phyetat0 ("startfi.nc",0,0,
345	& nsoilmx,nq,
346	& day_ini,time_phys,
347	& tsurf,tsoil,emis,q2,qsurf,co2ice)
348	#else
349	#include "meso_inc/meso_inc_ini.F"
350	#endif
351
352	if (pday.ne.day_ini) then
353	write(,) "PHYSIQ: ERROR: bad synchronization between ",
354	& "physics and dynamics"
355	write(,) "dynamics day: ",pday
356	write(,) "physics day: ",day_ini
357	stop
358	endif
359
360	write (,) 'In physiq day_ini =', day_ini
361
362	c Initialize albedo and orbital calculation
363	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
364	CALL surfini(ngrid,co2ice,qsurf,albedo)
365	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
366
367	c initialize soil
368	c ~~~~~~~~~~~~~~~
369	IF (callsoil) THEN
370	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
371	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
372	ELSE
373	PRINT*,
374	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
375	DO ig=1,ngrid
376	capcal(ig)=1.e5
377	fluxgrd(ig)=0.
378	ENDDO
379	ENDIF
380	icount=1
381
382
383	c initialize tracers
384	c ~~~~~~~~~~~~~~~~~~
385	tracerdyn=tracer
386	IF (tracer) THEN
387	CALL initracer(qsurf,co2ice)
388	ENDIF ! end tracer
389
390	#ifdef MESOSCALE
391	#include "meso_inc/meso_inc_caps.F"
392	#endif
393
394	#ifndef MESOSCALE
395	c Determining gridpoint near Viking Lander 1 (used for diagnostic only)
396	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
397
398	if(ngrid.ne.1) then
399	latvl1= 22.27
400	lonvl1= -47.94
401	ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim +
402	& int(1.5+(lonvl1+180)*iim/360.)
403	write(,) 'Viking Lander 1 GCM point: lat,lon',
404	& lati(ig_vl1)180/pi, long(ig_vl1)180/pi
405	end if
406	#endif
407
408	#ifndef MESOSCALE
409	c Initialize thermospheric parameters
410	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
411
412	if (callthermos) call param_read
413	#endif
414	c Initialize R and Cp as constant
415
416	if (.not.callthermos .and. .not.photochem) then
417	do l=1,nlayermx
418	do ig=1,ngridmx
419	rnew(ig,l)=r
420	cpnew(ig,l)=cpp
421	mmean(ig,l)=mugaz
422	enddo
423	enddo
424	endif
425
426	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
427	write(,)"physiq: water_param Surface ice alb:",alb_surfice
428	ENDIF
429
430	ENDIF ! (end of "if firstcall")
431
432	c ---------------------------------------------------
433	c 1.2 Initializations done at every physical timestep:
434	c ---------------------------------------------------
435	c
436	IF (ngrid.NE.ngridmx) THEN
437	PRINT*,'STOP in PHYSIQ'
438	PRINT*,'Probleme de dimensions :'
439	PRINT*,'ngrid = ',ngrid
440	PRINT*,'ngridmx = ',ngridmx
441	STOP
442	ENDIF
443
444	c Initialize various variables
445	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
446	call zerophys(ngrid*nlayer, pdv)
447	call zerophys(ngrid*nlayer, pdu)
448	call zerophys(ngrid*nlayer, pdt)
449	call zerophys(ngridnlayernq, pdq)
450	call zerophys(ngrid, pdpsrf)
451	call zerophys(ngrid, zflubid)
452	call zerophys(ngrid, zdtsurf)
453	call zerophys(ngrid*nq, dqsurf)
454	igout=ngrid/2+1
455
456
457	zday=pday+ptime ! compute time, in sols (and fraction thereof)
458
459	c Compute Solar Longitude (Ls) :
460	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
461	if (season) then
462	call solarlong(zday,zls)
463	else
464	call solarlong(float(day_ini),zls)
465	end if
466
467	c Compute geopotential at interlayers
468	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
469	c ponderation des altitudes au niveau des couches en dp/p
470
471	DO l=1,nlayer
472	DO ig=1,ngrid
473	zzlay(ig,l)=pphi(ig,l)/g
474	ENDDO
475	ENDDO
476	DO ig=1,ngrid
477	zzlev(ig,1)=0.
478	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
479	ENDDO
480	DO l=2,nlayer
481	DO ig=1,ngrid
482	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
483	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
484	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
485	ENDDO
486	ENDDO
487
488
489	! Potential temperature calculation not the same in physiq and dynamic
490
491	c Compute potential temperature
492	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
493	DO l=1,nlayer
494	DO ig=1,ngrid
495	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
496	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
497	ENDDO
498	ENDDO
499
500	#ifndef MESOSCALE
501	c-----------------------------------------------------------------------
502	c 1.2.5 Compute mean mass, cp, and R
503	c --------------------------------
504
505	if(photochem.or.callthermos) then
506	call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
507	endif
508	#endif
509	c-----------------------------------------------------------------------
510	c 2. Compute radiative tendencies :
511	c------------------------------------
512
513
514	IF (callrad) THEN
515	IF( MOD(icount-1,iradia).EQ.0) THEN
516
517	#ifdef MESOSCALE
518	write (,) 'call radiative transfer'
519	#endif
520	c Local Solar zenith angle
521	c ~~~~~~~~~~~~~~~~~~~~~~~~
522	CALL orbite(zls,dist_sol,declin)
523
524	IF(diurnal) THEN
525	ztim1=SIN(declin)
526	ztim2=COS(declin)COS(2.pi*(zday-.5))
527	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
528
529	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
530	s ztim1,ztim2,ztim3, mu0,fract)
531
532	ELSE
533	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
534	ENDIF
535
536	c NLTE cooling from CO2 emission
537	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
538
539	IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
540
541	c Find number of layers for LTE radiation calculations
542	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
543	& CALL nlthermeq(ngrid,nlayer,pplev,pplay)
544
545	c Note: Dustopacity.F has been transferred to callradite.F
546
547	c Call main radiative transfer scheme
548	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
549	c Transfer through CO2 (except NIR CO2 absorption)
550	c and aerosols (dust and water ice)
551
552	c Radiative transfer
553	c ------------------
554
555	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
556	$ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
557	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
558	& tauref,tau,aerosol,ccn,rdust,rice,nuice)
559
560	#ifdef MESOSCALE
561	#include "meso_inc/meso_inc_slope.F"
562	#endif
563
564	c CO2 near infrared absorption
565	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
566	call zerophys(ngrid*nlayer,zdtnirco2)
567	if (callnirco2) then
568	call nirco2abs (ngrid,nlayer,pplay,dist_sol,
569	. mu0,fract,declin, zdtnirco2)
570	endif
571
572	c Radiative flux from the sky absorbed by the surface (W.m-2)
573	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
574	DO ig=1,ngrid
575	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
576	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
577	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
578	ENDDO
579
580
581	c Net atmospheric radiative heating rate (K.s-1)
582	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
583	IF(callnlte) THEN
584	CALL blendrad(ngrid, nlayer, pplay,
585	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
586	ELSE
587	DO l=1,nlayer
588	DO ig=1,ngrid
589	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
590	& +zdtnirco2(ig,l)
591	ENDDO
592	ENDDO
593	ENDIF
594
595	ENDIF ! of if(mod(icount-1,iradia).eq.0)
596
597	c Transformation of the radiative tendencies:
598	c -------------------------------------------
599
600	c Net radiative surface flux (W.m-2)
601	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
602	c
603	DO ig=1,ngrid
604	zplanck(ig)=tsurf(ig)*tsurf(ig)
605	zplanck(ig)=emis(ig)*
606	$ stephanzplanck(ig)zplanck(ig)
607	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
608	#ifdef MESOSCALE
609	!!!! param slope
610	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
611	fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
612	#endif
613	ENDDO
614
615	DO l=1,nlayer
616	DO ig=1,ngrid
617	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
618	ENDDO
619	ENDDO
620
621	ENDIF ! of IF (callrad)
622
623	#ifndef MESOSCALE
624	c-----------------------------------------------------------------------
625	c 3. Gravity wave and subgrid scale topography drag :
626	c -------------------------------------------------
627
628
629	IF(calllott)THEN
630
631	CALL calldrag_noro(ngrid,nlayer,ptimestep,
632	& pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
633
634	DO l=1,nlayer
635	DO ig=1,ngrid
636	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
637	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
638	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
639	ENDDO
640	ENDDO
641	ENDIF
642	#endif
643	c-----------------------------------------------------------------------
644	c 4. Vertical diffusion (turbulent mixing):
645	c -----------------------------------------
646
647	IF (calldifv) THEN
648
649	DO ig=1,ngrid
650	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
651	ENDDO
652
653	CALL zerophys(ngrid*nlayer,zdum1)
654	CALL zerophys(ngrid*nlayer,zdum2)
655	do l=1,nlayer
656	do ig=1,ngrid
657	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
658	enddo
659	enddo
660
661	c Calling vdif (Martian version WITH CO2 condensation)
662	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
663	$ ptimestep,capcal,lwrite,
664	$ pplay,pplev,zzlay,zzlev,z0,
665	$ pu,pv,zh,pq,tsurf,emis,qsurf,
666	$ zdum1,zdum2,zdh,pdq,zflubid,
667	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
668	& zdqdif,zdqsdif)
669
670	#ifdef MESOSCALE
671	#include "meso_inc/meso_inc_les.F"
672	#endif
673	DO l=1,nlayer
674	DO ig=1,ngrid
675	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
676	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
677	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
678
679	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
680
681	ENDDO
682	ENDDO
683
684	DO ig=1,ngrid
685	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
686	ENDDO
687
688	if (tracer) then
689	DO iq=1, nq
690	DO l=1,nlayer
691	DO ig=1,ngrid
692	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
693	ENDDO
694	ENDDO
695	ENDDO
696	DO iq=1, nq
697	DO ig=1,ngrid
698	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
699	ENDDO
700	ENDDO
701	end if ! of if (tracer)
702
703	ELSE
704	DO ig=1,ngrid
705	zdtsurf(ig)=zdtsurf(ig)+
706	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
707	ENDDO
708	#ifdef MESOSCALE
709	IF (flag_LES) THEN
710	write(,) 'LES mode !'
711	write(,) 'Please set calldifv to T in callphys.def'
712	STOP
713	ENDIF
714	#endif
715	ENDIF ! of IF (calldifv)
716
717	c-----------------------------------------------------------------------
718	c TEST. Thermals :
719	c HIGHLY EXPERIMENTAL, BEWARE !!
720	c -----------------------------
721
722	if(calltherm) then
723
724	call calltherm_interface(firstcall,
725	$ long,lati,zzlev,zzlay,
726	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
727	$ pplay,pplev,pphi,zpopsk,
728	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
729	$ dtke_th,hfmax_th,wmax_th)
730
731	DO l=1,nlayer
732	DO ig=1,ngrid
733	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
734	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
735	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
736	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
737	ENDDO
738	ENDDO
739
740	DO ig=1,ngrid
741	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
742	ENDDO
743
744	if (tracer) then
745	DO iq=1,nq
746	DO l=1,nlayer
747	DO ig=1,ngrid
748	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
749	ENDDO
750	ENDDO
751	ENDDO
752	endif
753
754	else !of if calltherm
755	lmax_th(:)=0
756	end if
757
758	c-----------------------------------------------------------------------
759	c 5. Dry convective adjustment:
760	c -----------------------------
761
762	IF(calladj) THEN
763
764	DO l=1,nlayer
765	DO ig=1,ngrid
766	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
767	ENDDO
768	ENDDO
769	CALL zerophys(ngrid*nlayer,zduadj)
770	CALL zerophys(ngrid*nlayer,zdvadj)
771	CALL zerophys(ngrid*nlayer,zdhadj)
772
773	CALL convadj(ngrid,nlayer,nq,ptimestep,
774	$ pplay,pplev,zpopsk,lmax_th,
775	$ pu,pv,zh,pq,
776	$ pdu,pdv,zdh,pdq,
777	$ zduadj,zdvadj,zdhadj,
778	$ zdqadj)
779
780
781	DO l=1,nlayer
782	DO ig=1,ngrid
783	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
784	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
785	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
786
787	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
788	ENDDO
789	ENDDO
790
791	if(tracer) then
792	DO iq=1, nq
793	DO l=1,nlayer
794	DO ig=1,ngrid
795	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
796	ENDDO
797	ENDDO
798	ENDDO
799	end if
800	ENDIF ! of IF(calladj)
801
802	c-----------------------------------------------------------------------
803	c 6. Carbon dioxide condensation-sublimation:
804	c -------------------------------------------
805
806	IF (callcond) THEN
807	CALL newcondens(ngrid,nlayer,nq,ptimestep,
808	$ capcal,pplay,pplev,tsurf,pt,
809	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
810	$ co2ice,albedo,emis,
811	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
812	$ fluxsurf_sw,zls)
813
814	DO l=1,nlayer
815	DO ig=1,ngrid
816	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
817	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
818	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
819	ENDDO
820	ENDDO
821	DO ig=1,ngrid
822	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
823	ENDDO
824
825	IF (tracer) THEN
826	DO iq=1, nq
827	DO l=1,nlayer
828	DO ig=1,ngrid
829	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
830	ENDDO
831	ENDDO
832	ENDDO
833	ENDIF ! of IF (tracer)
834
835	ENDIF ! of IF (callcond)
836
837	c-----------------------------------------------------------------------
838	c 7. Specific parameterizations for tracers
839	c: -----------------------------------------
840
841	if (tracer) then
842
843	c 7a. Water and ice
844	c ---------------
845
846	c ---------------------------------------
847	c Water ice condensation in the atmosphere
848	c ----------------------------------------
849	IF (water) THEN
850
851	call watercloud(ngrid,nlayer,ptimestep,
852	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
853	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
854	& nq,naerkind,tau,
855	& ccn,rdust,rice,nuice)
856	if (activice) then
857	c Temperature variation due to latent heat release
858	DO l=1,nlayer
859	DO ig=1,ngrid
860	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
861	ENDDO
862	ENDDO
863	endif
864
865	! increment water vapour and ice atmospheric tracers tendencies
866	IF (water) THEN
867	DO l=1,nlayer
868	DO ig=1,ngrid
869	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
870	& zdqcloud(ig,l,igcm_h2o_vap)
871	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
872	& zdqcloud(ig,l,igcm_h2o_ice)
873	ENDDO
874	ENDDO
875	ENDIF ! of IF (water) THEN
876	! Increment water ice surface tracer tendency
877	DO ig=1,ngrid
878	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
879	& zdqscloud(ig,igcm_h2o_ice)
880	ENDDO
881
882	END IF ! of IF (water)
883
884	c 7b. Chemical species
885	c ------------------
886
887	#ifndef MESOSCALE
888	c --------------
889	c photochemistry :
890	c --------------
891	IF (photochem .or. thermochem) then
892	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
893	& zzlay,zday,pq,pdq,rice,
894	& zdqchim,zdqschim,zdqcloud,zdqscloud)
895	!NB: Photochemistry includes condensation of H2O2
896
897	! increment values of tracers:
898	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
899	! tracers is zero anyways
900	DO l=1,nlayer
901	DO ig=1,ngrid
902	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
903	ENDDO
904	ENDDO
905	ENDDO ! of DO iq=1,nq
906	! add condensation tendency for H2O2
907	if (igcm_h2o2.ne.0) then
908	DO l=1,nlayer
909	DO ig=1,ngrid
910	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
911	& +zdqcloud(ig,l,igcm_h2o2)
912	ENDDO
913	ENDDO
914	endif
915
916	! increment surface values of tracers:
917	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
918	! tracers is zero anyways
919	DO ig=1,ngrid
920	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
921	ENDDO
922	ENDDO ! of DO iq=1,nq
923	! add condensation tendency for H2O2
924	if (igcm_h2o2.ne.0) then
925	DO ig=1,ngrid
926	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
927	& +zdqscloud(ig,igcm_h2o2)
928	ENDDO
929	endif
930
931	END IF ! of IF (photochem.or.thermochem)
932	#endif
933
934	c 7c. Aerosol particles
935	c -------------------
936
937	c ----------
938	c Dust devil :
939	c ----------
940	IF(callddevil) then
941	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
942	& zdqdev,zdqsdev)
943
944	if (dustbin.ge.1) then
945	do iq=1,nq
946	DO l=1,nlayer
947	DO ig=1,ngrid
948	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
949	ENDDO
950	ENDDO
951	enddo
952	do iq=1,nq
953	DO ig=1,ngrid
954	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
955	ENDDO
956	enddo
957	endif ! of if (dustbin.ge.1)
958
959	END IF ! of IF (callddevil)
960
961	c -------------
962	c Sedimentation : acts also on water ice
963	c -------------
964	IF (sedimentation) THEN
965	!call zerophys(ngridnlayernq, zdqsed)
966	zdqsed(1:ngrid,1:nlayer,1:nq)=0
967	!call zerophys(ngrid*nq, zdqssed)
968	zdqssed(1:ngrid,1:nq)=0
969
970	call callsedim(ngrid,nlayer, ptimestep,
971	& pplev,zzlev, pt, rdust, rice,
972	& pq, pdq, zdqsed, zdqssed,nq)
973
974	DO iq=1, nq
975	DO l=1,nlayer
976	DO ig=1,ngrid
977	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
978	ENDDO
979	ENDDO
980	ENDDO
981	DO iq=1, nq
982	DO ig=1,ngrid
983	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
984	ENDDO
985	ENDDO
986	END IF ! of IF (sedimentation)
987
988	c 7d. Updates
989	c ---------
990
991	DO iq=1, nq
992	DO ig=1,ngrid
993
994	c ---------------------------------
995	c Updating tracer budget on surface
996	c ---------------------------------
997	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
998
999	ENDDO ! (ig)
1000	ENDDO ! (iq)
1001
1002	endif ! of if (tracer)
1003
1004	#ifndef MESOSCALE
1005	c-----------------------------------------------------------------------
1006	c 8. THERMOSPHERE CALCULATION
1007	c-----------------------------------------------------------------------
1008
1009	if (callthermos) then
1010	call thermosphere(pplev,pplay,dist_sol,
1011	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1012	& pt,pq,pu,pv,pdt,pdq,
1013	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1014
1015	DO l=1,nlayer
1016	DO ig=1,ngrid
1017	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1018	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1019	& +zdteuv(ig,l)
1020	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1021	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1022	DO iq=1, nq
1023	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1024	ENDDO
1025	ENDDO
1026	ENDDO
1027
1028	endif ! of if (callthermos)
1029	#endif
1030
1031	c-----------------------------------------------------------------------
1032	c 9. Surface and sub-surface soil temperature
1033	c-----------------------------------------------------------------------
1034	c
1035	c
1036	c 9.1 Increment Surface temperature:
1037	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1038
1039	DO ig=1,ngrid
1040	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1041	ENDDO
1042
1043	c Prescribe a cold trap at south pole (except at high obliquity !!)
1044	c Temperature at the surface is set there to be the temperature
1045	c corresponding to equilibrium temperature between phases of CO2
1046
1047	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
1048	#ifndef MESOSCALE
1049	if (caps.and.(obliquit.lt.27.)) then
1050	! NB: Updated surface pressure, at grid point 'ngrid', is
1051	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1052	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1053	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1054	endif
1055	#endif
1056	c -------------------------------------------------------------
1057	c Change of surface albedo (set to 0.4) in case of ground frost
1058	c everywhere except on the north permanent cap and in regions
1059	c covered by dry ice.
1060	c ALWAYS PLACE these lines after newcondens !!!
1061	c -------------------------------------------------------------
1062	do ig=1,ngrid
1063	if ((co2ice(ig).eq.0).and.
1064	& (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
1065	albedo(ig,1) = alb_surfice
1066	albedo(ig,2) = alb_surfice
1067	endif
1068	enddo ! of do ig=1,ngrid
1069	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
1070
1071	c
1072	c 9.2 Compute soil temperatures and subsurface heat flux:
1073	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1074	IF (callsoil) THEN
1075	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1076	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
1077	ENDIF
1078
1079	c-----------------------------------------------------------------------
1080	c 10. Write output files
1081	c ----------------------
1082
1083	c -------------------------------
1084	c Dynamical fields incrementation
1085	c -------------------------------
1086	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1087	! temperature, zonal and meridional wind
1088	DO l=1,nlayer
1089	DO ig=1,ngrid
1090	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1091	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1092	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1093	ENDDO
1094	ENDDO
1095
1096	! tracers
1097	DO iq=1, nq
1098	DO l=1,nlayer
1099	DO ig=1,ngrid
1100	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1101	ENDDO
1102	ENDDO
1103	ENDDO
1104
1105	! surface pressure
1106	DO ig=1,ngrid
1107	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1108	ENDDO
1109
1110	! pressure
1111	DO l=1,nlayer
1112	DO ig=1,ngrid
1113	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1114	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1115	ENDDO
1116	ENDDO
1117
1118	! Density
1119	DO l=1,nlayer
1120	DO ig=1,ngrid
1121	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1122	ENDDO
1123	ENDDO
1124
1125	c Compute surface stress : (NB: z0 is a common in surfdat.h)
1126	c DO ig=1,ngrid
1127	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
1128	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1129	c ENDDO
1130
1131	c Sum of fluxes in solar spectral bands (for output only)
1132	DO ig=1,ngrid
1133	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1134	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1135	ENDDO
1136	c ***** TEST ****************************************************
1137	ztim1 = 999
1138	DO l=1,nlayer
1139	DO ig=1,ngrid
1140	if (pt(ig,l).lt.ztim1) then
1141	ztim1 = pt(ig,l)
1142	igmin = ig
1143	lmin = l
1144	end if
1145	ENDDO
1146	ENDDO
1147	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1148	write(,) 'PHYSIQ: stability WARNING :'
1149	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1150	& 'ig l =', igmin, lmin
1151	end if
1152	c *******************************************************************
1153
1154	c ---------------------
1155	c Outputs to the screen
1156	c ---------------------
1157
1158	IF (lwrite) THEN
1159	PRINT*,'Global diagnostics for the physics'
1160	PRINT,'Variables and their increments x and dx/dt dt'
1161	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1162	WRITE(*,'(2f10.5,2f15.5)')
1163	s tsurf(igout),zdtsurf(igout)*ptimestep,
1164	s pplev(igout,1),pdpsrf(igout)*ptimestep
1165	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1166	WRITE(*,'(i4,6f10.5)') (l,
1167	s pu(igout,l),pdu(igout,l)*ptimestep,
1168	s pv(igout,l),pdv(igout,l)*ptimestep,
1169	s pt(igout,l),pdt(igout,l)*ptimestep,
1170	s l=1,nlayer)
1171	ENDIF ! of IF (lwrite)
1172
1173	IF (ngrid.NE.1) THEN
1174	print,'Ls =',zls180./pi
1175	& ,' tauref(700 Pa,lat=0) =',tauref(ngrid/2)
1176	#ifndef MESOSCALE
1177	& ,' tau(Viking1) =',tau(ig_vl1,1)
1178	#endif
1179
1180	#ifndef MESOSCALE
1181	c -------------------------------------------------------------------
1182	c Writing NetCDF file "RESTARTFI" at the end of the run
1183	c -------------------------------------------------------------------
1184	c Note: 'restartfi' is stored just before dynamics are stored
1185	c in 'restart'. Between now and the writting of 'restart',
1186	c there will have been the itau=itau+1 instruction and
1187	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1188	c thus we store for time=time+dtvr
1189
1190	IF(lastcall) THEN
1191	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1192	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1193	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1194	. ptimestep,pday,
1195	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
1196	. area,albedodat,inertiedat,zmea,zstd,zsig,
1197	. zgam,zthe)
1198	ENDIF
1199	#endif
1200
1201	c -------------------------------------------------------------------
1202	c Calculation of diagnostic variables written in both stats and
1203	c diagfi files
1204	c -------------------------------------------------------------------
1205
1206	if (tracer) then
1207	if (water) then
1208
1209	call zerophys(ngrid,mtot)
1210	call zerophys(ngrid,icetot)
1211	call zerophys(ngrid,rave)
1212	call zerophys(ngrid,tauTES)
1213	do ig=1,ngrid
1214	do l=1,nlayermx
1215	mtot(ig) = mtot(ig) +
1216	& zq(ig,l,igcm_h2o_vap) *
1217	& (pplev(ig,l) - pplev(ig,l+1)) / g
1218	icetot(ig) = icetot(ig) +
1219	& zq(ig,l,igcm_h2o_ice) *
1220	& (pplev(ig,l) - pplev(ig,l+1)) / g
1221	rave(ig) = rave(ig) +
1222	& zq(ig,l,igcm_h2o_ice) *
1223	& (pplev(ig,l) - pplev(ig,l+1)) / g *
1224	& rice(ig,l) * (1.+nuice_ref)
1225	c Computing abs optical depth at 825 cm-1 in each
1226	c layer to simulate NEW TES retrieval
1227	Qabsice = min(
1228	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
1229	& )
1230	opTES(ig,l)= 0.75 * Qabsice *
1231	& zq(ig,l,igcm_h2o_ice) *
1232	& (pplev(ig,l) - pplev(ig,l+1)) / g
1233	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1234	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1235	enddo
1236	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
1237	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1238	enddo
1239
1240	endif ! of if (water)
1241	endif ! of if (tracer)
1242
1243	c -----------------------------------------------------------------
1244	c WSTATS: Saving statistics
1245	c -----------------------------------------------------------------
1246	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1247	c which can later be used to make the statistic files of the run:
1248	c "stats") only possible in 3D runs !
1249
1250	#ifndef MESOSCALE
1251	IF (callstats) THEN
1252
1253	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1254	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1255	call wstats(ngrid,"co2ice","CO2 ice cover",
1256	& "kg.m-2",2,co2ice)
1257	call wstats(ngrid,"fluxsurf_lw",
1258	& "Thermal IR radiative flux to surface","W.m-2",2,
1259	& fluxsurf_lw)
1260	call wstats(ngrid,"fluxsurf_sw",
1261	& "Solar radiative flux to surface","W.m-2",2,
1262	& fluxsurf_sw_tot)
1263	call wstats(ngrid,"fluxtop_lw",
1264	& "Thermal IR radiative flux to space","W.m-2",2,
1265	& fluxtop_lw)
1266	call wstats(ngrid,"fluxtop_sw",
1267	& "Solar radiative flux to space","W.m-2",2,
1268	& fluxtop_sw_tot)
1269	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1270	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1271	call wstats(ngrid,"v","Meridional (North-South) wind",
1272	& "m.s-1",3,zv)
1273	call wstats(ngrid,"w","Vertical (down-up) wind",
1274	& "m.s-1",3,pw)
1275	call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
1276	call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
1277	c call wstats(ngrid,"q2",
1278	c & "Boundary layer eddy kinetic energy",
1279	c & "m2.s-2",3,q2)
1280	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1281	c & emis)
1282	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
1283	c & 2,zstress)
1284	c call wstats(ngrid,"sw_htrt","sw heat.rate",
1285	c & "W.m-2",3,zdtsw)
1286	c call wstats(ngrid,"lw_htrt","lw heat.rate",
1287	c & "W.m-2",3,zdtlw)
1288
1289	if (tracer) then
1290	if (water) then
1291	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1292	& *mugaz/mmol(igcm_h2o_vap)
1293	call wstats(ngrid,"vmr_h2ovapor",
1294	& "H2O vapor volume mixing ratio","mol/mol",
1295	& 3,vmr)
1296	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1297	& *mugaz/mmol(igcm_h2o_ice)
1298	call wstats(ngrid,"vmr_h2oice",
1299	& "H2O ice volume mixing ratio","mol/mol",
1300	& 3,vmr)
1301	call wstats(ngrid,"h2o_ice_s",
1302	& "surface h2o_ice","kg/m2",
1303	& 2,qsurf(1,igcm_h2o_ice))
1304
1305	call wstats(ngrid,"mtot",
1306	& "total mass of water vapor","kg/m2",
1307	& 2,mtot)
1308	call wstats(ngrid,"icetot",
1309	& "total mass of water ice","kg/m2",
1310	& 2,icetot)
1311	call wstats(ngrid,"reffice",
1312	& "Mean reff","m",
1313	& 2,rave)
1314	c call wstats(ngrid,"rice",
1315	c & "Ice particle size","m",
1316	c & 3,rice)
1317	c If activice is true, tauTES is computed in aeropacity.F;
1318	if (.not.activice) then
1319	call wstats(ngrid,"tauTESap",
1320	& "tau abs 825 cm-1","",
1321	& 2,tauTES)
1322	endif
1323
1324	endif ! of if (water)
1325
1326	if (thermochem.or.photochem) then
1327	do iq=1,nq
1328	if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
1329	. (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
1330	. (noms(iq).eq."h2").or.
1331	. (noms(iq).eq."o3")) then
1332	do l=1,nlayer
1333	do ig=1,ngrid
1334	vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
1335	end do
1336	end do
1337	call wstats(ngrid,"vmr_"//trim(noms(iq)),
1338	. "Volume mixing ratio","mol/mol",3,vmr)
1339	endif
1340	enddo
1341	endif ! of if (thermochem.or.photochem)
1342
1343	endif ! of if (tracer)
1344
1345	IF(lastcall) THEN
1346	write (,) "Writing stats..."
1347	call mkstats(ierr)
1348	ENDIF
1349
1350	ENDIF !if callstats
1351	#endif
1352
1353	c (Store EOF for Mars Climate database software)
1354	IF (calleofdump) THEN
1355	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1356	ENDIF
1357
1358	#ifdef MESOSCALE
1359	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
1360	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
1361	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
1362	IF (igcm_h2o_ice .ne. 0) THEN
1363	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
1364	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1365	. * mugaz / mmol(igcm_h2o_ice)
1366	ENDIF
1367	c AUTOMATICALLY GENERATED FROM REGISTRY
1368	#include "fill_save.inc"
1369	#else
1370
1371	c ==========================================================
1372	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
1373	c any variable for diagnostic (output with period
1374	c "ecritphy", set in "run.def")
1375	c ==========================================================
1376	c WRITEDIAGFI can ALSO be called from any other subroutines
1377	c for any variables !!
1378	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1379	c & emis)
1380	! call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
1381	! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
1382	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1383	& tsurf)
1384	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1385	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
1386	& co2ice)
1387	c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
1388	c & "K",2,zt(1,7))
1389	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1390	& fluxsurf_lw)
1391	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1392	& fluxsurf_sw_tot)
1393	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1394	& fluxtop_lw)
1395	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1396	& fluxtop_sw_tot)
1397	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1398	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1399	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1400	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1401	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1402	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1403	! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1404	c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
1405	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
1406	c & zstress)
1407	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
1408	c & 'w.m-2',3,zdtsw)
1409	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
1410	c & 'w.m-2',3,zdtlw)
1411	c CALL WRITEDIAGFI(ngridmx,'tauTESap',
1412	c & 'tau abs 825 cm-1',
1413	c & '',2,tauTES)
1414
1415
1416	c ----------------------------------------------------------
1417	c Outputs of the CO2 cycle
1418	c ----------------------------------------------------------
1419
1420	if (tracer.and.(igcm_co2.ne.0)) then
1421	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
1422	! & "kg/kg",2,zq(1,1,igcm_co2))
1423	! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
1424	! & "kg/kg",3,zq(1,1,igcm_co2))
1425
1426	! Compute co2 column
1427	call zerophys(ngrid,co2col)
1428	do l=1,nlayermx
1429	do ig=1,ngrid
1430	co2col(ig)=co2col(ig)+
1431	& zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
1432	enddo
1433	enddo
1434	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
1435	& co2col)
1436	endif ! of if (tracer.and.(igcm_co2.ne.0))
1437
1438	c ----------------------------------------------------------
1439	c Outputs of the water cycle
1440	c ----------------------------------------------------------
1441	if (tracer) then
1442	if (water) then
1443
1444	CALL WRITEDIAGFI(ngridmx,'mtot',
1445	& 'total mass of water vapor',
1446	& 'kg/m2',2,mtot)
1447	CALL WRITEDIAGFI(ngridmx,'icetot',
1448	& 'total mass of water ice',
1449	& 'kg/m2',2,icetot)
1450	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1451	c & *mugaz/mmol(igcm_h2o_ice)
1452	c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
1453	c & 'mol/mol',3,vmr)
1454	CALL WRITEDIAGFI(ngridmx,'reffice',
1455	& 'Mean reff',
1456	& 'm',2,rave)
1457	c call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
1458	c & 'm',3,rice)
1459	c If activice is true, tauTES is computed in aeropacity.F;
1460	if (.not.activice) then
1461	CALL WRITEDIAGFI(ngridmx,'tauTESap',
1462	& 'tau abs 825 cm-1',
1463	& '',2,tauTES)
1464	endif
1465	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
1466	& 'surface h2o_ice',
1467	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
1468	endif !(water)
1469
1470
1471	if (water.and..not.photochem) then
1472	iq=nq
1473	c write(str2(1:2),'(i2.2)') iq
1474	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
1475	c & 'kg.m-2',2,zdqscloud(1,iq))
1476	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
1477	c & 'kg/kg',3,zdqchim(1,1,iq))
1478	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
1479	c & 'kg/kg',3,zdqdif(1,1,iq))
1480	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
1481	c & 'kg/kg',3,zdqadj(1,1,iq))
1482	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
1483	c & 'kg/kg',3,zdqc(1,1,iq))
1484	endif !(water.and..not.photochem)
1485	endif
1486
1487	c ----------------------------------------------------------
1488	c Outputs of the dust cycle
1489	c ----------------------------------------------------------
1490
1491	c call WRITEDIAGFI(ngridmx,'tauref',
1492	c & 'Dust ref opt depth','NU',2,tauref)
1493
1494	if (tracer.and.(dustbin.ne.0)) then
1495	c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
1496	if (doubleq) then
1497	c call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
1498	c & 'kg.m-2',2,qsurf(1,1))
1499	c call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
1500	c & 'N.m-2',2,qsurf(1,2))
1501	c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
1502	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
1503	c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
1504	c & 'kg.m-2.s-1',2,zdqssed(1,1))
1505	call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
1506	& 'm',3,rdust*ref_r0)
1507	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
1508	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
1509	c call WRITEDIAGFI(ngridmx,'dustN','Dust number',
1510	c & 'part/kg',3,pq(1,1,igcm_dust_number))
1511	else
1512	do iq=1,dustbin
1513	write(str2(1:2),'(i2.2)') iq
1514	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
1515	& 'kg/kg',3,zq(1,1,iq))
1516	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
1517	& 'kg.m-2',2,qsurf(1,iq))
1518	end do
1519	endif ! (doubleq)
1520	c if (submicron) then
1521	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
1522	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
1523	c endif ! (submicron)
1524	end if ! (tracer.and.(dustbin.ne.0))
1525
1526	c ----------------------------------------------------------
1527	c Outputs of thermals
1528	c ----------------------------------------------------------
1529	if (calltherm) then
1530
1531	! call WRITEDIAGFI(ngrid,'dtke',
1532	! & 'tendance tke thermiques','m2/s2',
1533	! & 3,dtke_th)
1534	! call WRITEDIAGFI(ngrid,'d_u_ajs',
1535	! & 'tendance u thermiques','m/s',
1536	! & 3,pdu_th*ptimestep)
1537	! call WRITEDIAGFI(ngrid,'d_v_ajs',
1538	! & 'tendance v thermiques','m/s',
1539	! & 3,pdv_th*ptimestep)
1540	! if (tracer) then
1541	! if (nq .eq. 2) then
1542	! call WRITEDIAGFI(ngrid,'deltaq_th',
1543	! & 'delta q thermiques','kg/kg',
1544	! & 3,ptimestep*pdq_th(:,:,2))
1545	! endif
1546	! endif
1547
1548	lmax_th_out(:)=real(lmax_th(:))
1549
1550	call WRITEDIAGFI(ngridmx,'lmax_th',
1551	& 'hauteur du thermique','K',
1552	& 2,lmax_th_out)
1553	call WRITEDIAGFI(ngridmx,'hfmax_th',
1554	& 'maximum TH heat flux','K.m/s',
1555	& 2,hfmax_th)
1556	call WRITEDIAGFI(ngridmx,'wmax_th',
1557	& 'maximum TH vertical velocity','m/s',
1558	& 2,wmax_th)
1559
1560	endif
1561
1562	c ----------------------------------------------------------
1563	c Output in netcdf file "diagsoil.nc" for subterranean
1564	c variables (output every "ecritphy", as for writediagfi)
1565	c ----------------------------------------------------------
1566
1567	! Write soil temperature
1568	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
1569	! & 3,tsoil)
1570	! Write surface temperature
1571	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
1572	! & 2,tsurf)
1573
1574	c ==========================================================
1575	c END OF WRITEDIAGFI
1576	c ==========================================================
1577	#endif
1578
1579	ELSE ! if(ngrid.eq.1)
1580
1581	#ifndef MESOSCALE
1582	print,'Ls =',zls180./pi,
1583	& ' tauref(700 Pa) =',tauref
1584	c ----------------------------------------------------------------------
1585	c Output in grads file "g1d" (ONLY when using testphys1d)
1586	c (output at every X physical timestep)
1587	c ----------------------------------------------------------------------
1588	c
1589	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
1590	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
1591	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
1592
1593	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
1594	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
1595	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
1596	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
1597
1598	! THERMALS STUFF 1D
1599
1600	if(calltherm) then
1601
1602	lmax_th_out(:)=real(lmax_th(:))
1603
1604	call WRITEDIAGFI(ngridmx,'lmax_th',
1605	& 'hauteur du thermique','point',
1606	& 0,lmax_th_out)
1607	call WRITEDIAGFI(ngridmx,'hfmax_th',
1608	& 'maximum TH heat flux','K.m/s',
1609	& 0,hfmax_th)
1610	call WRITEDIAGFI(ngridmx,'wmax_th',
1611	& 'maximum TH vertical velocity','m/s',
1612	& 0,wmax_th)
1613
1614
1615	co2col(:)=0.
1616	dummycol(:)=0.
1617	if (tracer) then
1618	do l=1,nlayermx
1619	do ig=1,ngrid
1620	co2col(ig)=co2col(ig)+
1621	& zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
1622	if (nqmx .gt. 1) then
1623	iq=2 ! to avoid out-of-bounds spotting by picky compilers
1624	! when gcm is compiled with only one tracer
1625	dummycol(ig)=dummycol(ig)+
1626	& zq(ig,l,iq)*(pplev(ig,l)-pplev(ig,l+1))/g
1627	endif
1628	enddo
1629	enddo
1630
1631	end if
1632	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
1633	& ,'kg/m-2',0,co2col)
1634	call WRITEDIAGFI(ngrid,'dummycol','integrated dummy mass' &
1635	& ,'kg/m-2',0,dummycol)
1636	endif
1637	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
1638	& ,'m/s',1,pw)
1639	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
1640	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
1641	& tsurf)
1642
1643	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
1644	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
1645	! or output in diagfi.nc (for testphys1d)
1646	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
1647	call WRITEDIAGFI(ngridmx,'temp','Temperature',
1648	& 'K',1,zt)
1649
1650	if(tracer) then
1651	c CALL writeg1d(ngrid,1,tau,'tau','SI')
1652	do iq=1,nq
1653	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
1654	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
1655	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
1656	end do
1657	end if
1658
1659	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
1660
1661	do l=2,nlayer-1
1662	tmean=zt(1,l)
1663	if(zt(1,l).ne.zt(1,l-1))
1664	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
1665	zlocal(l)= zlocal(l-1)
1666	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
1667	enddo
1668	zlocal(nlayer)= zlocal(nlayer-1)-
1669	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
1670	& rnew(1,nlayer)*tmean/g
1671	#endif
1672
1673	END IF ! if(ngrid.ne.1)
1674
1675	icount=icount+1
1676	#ifdef MESOSCALE
1677	write(,) 'now, back to the dynamical core...'
1678	#endif
1679	RETURN
1680	END

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

Download in other formats:

Original Format