Context Navigation

physiq.F @ 528

Last change on this file since 528 was 528, checked in by aslmd, 14 years ago
LMDZ.MARS : problem with r520. commits were made on an older version. fix the (numerous) problems this caused. please developers use svn update before committing.
File size: 77.6 KB

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

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

Download in other formats:

Original Format