Context Navigation

physiq.F @ 629

Last change on this file since 629 was 627, checked in by emillour, 13 years ago

Mars GCM:

Some cleanup on messages and comments in code about the reference pressure

for dust opacity which is now 610 Pa.

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

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

Download in other formats:

Original Format