Context Navigation

physiq.F @ 623

Last change on this file since 623 was 608, checked in by aslmd, 13 years ago
LMDZ.MARS: minor correction in solarlong. MESOSCALE: output mtot and icetot now same units, notes updated, improved save_this_simu so that all .def are saved. UTIL PYTHON: now in mesoscale mode, sections show lat/lon instead of nx or ny (to be improved), also corrected 1D labels that was not working for several time asked.
File size: 80.1 KB

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

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

Download in other formats:

Original Format