Context Navigation

physiq.F @ 277

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

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

Download in other formats:

Original Format