Context Navigation

physiq.F @ 263

Last change on this file since 263 was 263, checked in by aslmd, 13 years ago
MESOSCALE: Something was broken for MESOSCALE model in r259-r260 (flag_LES is not known in vdifc). Moved this to meso_inc/meso_inc_les.F. Also there was a problem with the new surface layer if calltherm is false because then wmax has no value... Fixed.
File size: 60.3 KB

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

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

Download in other formats:

Original Format