Context Navigation

physiq.F @ 272

Last change on this file since 272 was 269, checked in by acolaitis, 14 years ago
Bug correction in call to surflayer, added recomputation of potential temperature zh after incrementation of physical tendancies.
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	else !of if calltherm
791	lmax_th(:)=0
792	wmax_th(:)=0.
793	end if
794
795	c-----------------------------------------------------------------------
796	c 5. Dry convective adjustment:
797	c -----------------------------
798
799	IF(calladj) THEN
800
801	DO l=1,nlayer
802	DO ig=1,ngrid
803	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
804	ENDDO
805	ENDDO
806	CALL zerophys(ngrid*nlayer,zduadj)
807	CALL zerophys(ngrid*nlayer,zdvadj)
808	CALL zerophys(ngrid*nlayer,zdhadj)
809
810	CALL convadj(ngrid,nlayer,nq,ptimestep,
811	$ pplay,pplev,zpopsk,lmax_th,
812	$ pu,pv,zh,pq,
813	$ pdu,pdv,zdh,pdq,
814	$ zduadj,zdvadj,zdhadj,
815	$ zdqadj)
816
817
818	DO l=1,nlayer
819	DO ig=1,ngrid
820	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
821	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
822	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
823
824	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
825	ENDDO
826	ENDDO
827
828	if(tracer) then
829	DO iq=1, nq
830	DO l=1,nlayer
831	DO ig=1,ngrid
832	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
833	ENDDO
834	ENDDO
835	ENDDO
836	end if
837	ENDIF ! of IF(calladj)
838
839	c-----------------------------------------------------------------------
840	c 6. Carbon dioxide condensation-sublimation:
841	c -------------------------------------------
842
843	IF (callcond) THEN
844	CALL newcondens(ngrid,nlayer,nq,ptimestep,
845	$ capcal,pplay,pplev,tsurf,pt,
846	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
847	$ co2ice,albedo,emis,
848	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
849	$ fluxsurf_sw,zls)
850
851	DO l=1,nlayer
852	DO ig=1,ngrid
853	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
854	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
855	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
856	ENDDO
857	ENDDO
858	DO ig=1,ngrid
859	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
860	ENDDO
861
862	IF (tracer) THEN
863	DO iq=1, nq
864	DO l=1,nlayer
865	DO ig=1,ngrid
866	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
867	ENDDO
868	ENDDO
869	ENDDO
870	ENDIF ! of IF (tracer)
871
872	ENDIF ! of IF (callcond)
873
874	c-----------------------------------------------------------------------
875	c 7. Specific parameterizations for tracers
876	c: -----------------------------------------
877
878	if (tracer) then
879
880	c 7a. Water and ice
881	c ---------------
882
883	c ---------------------------------------
884	c Water ice condensation in the atmosphere
885	c ----------------------------------------
886	IF (water) THEN
887
888	call watercloud(ngrid,nlayer,ptimestep,
889	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
890	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
891	& nq,naerkind,tau,
892	& ccn,rdust,rice,nuice)
893	if (activice) then
894	c Temperature variation due to latent heat release
895	DO l=1,nlayer
896	DO ig=1,ngrid
897	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
898	ENDDO
899	ENDDO
900	endif
901
902	! increment water vapour and ice atmospheric tracers tendencies
903	IF (water) THEN
904	DO l=1,nlayer
905	DO ig=1,ngrid
906	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
907	& zdqcloud(ig,l,igcm_h2o_vap)
908	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
909	& zdqcloud(ig,l,igcm_h2o_ice)
910	ENDDO
911	ENDDO
912	ENDIF ! of IF (water) THEN
913	! Increment water ice surface tracer tendency
914	DO ig=1,ngrid
915	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
916	& zdqscloud(ig,igcm_h2o_ice)
917	ENDDO
918
919	END IF ! of IF (water)
920
921	c 7b. Chemical species
922	c ------------------
923
924	#ifndef MESOSCALE
925	c --------------
926	c photochemistry :
927	c --------------
928	IF (photochem .or. thermochem) then
929	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
930	& zzlay,zday,pq,pdq,rice,
931	& zdqchim,zdqschim,zdqcloud,zdqscloud)
932	!NB: Photochemistry includes condensation of H2O2
933
934	! increment values of tracers:
935	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
936	! tracers is zero anyways
937	DO l=1,nlayer
938	DO ig=1,ngrid
939	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
940	ENDDO
941	ENDDO
942	ENDDO ! of DO iq=1,nq
943	! add condensation tendency for H2O2
944	if (igcm_h2o2.ne.0) then
945	DO l=1,nlayer
946	DO ig=1,ngrid
947	pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
948	& +zdqcloud(ig,l,igcm_h2o2)
949	ENDDO
950	ENDDO
951	endif
952
953	! increment surface values of tracers:
954	DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
955	! tracers is zero anyways
956	DO ig=1,ngrid
957	dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
958	ENDDO
959	ENDDO ! of DO iq=1,nq
960	! add condensation tendency for H2O2
961	if (igcm_h2o2.ne.0) then
962	DO ig=1,ngrid
963	dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
964	& +zdqscloud(ig,igcm_h2o2)
965	ENDDO
966	endif
967
968	END IF ! of IF (photochem.or.thermochem)
969	#endif
970
971	c 7c. Aerosol particles
972	c -------------------
973
974	c ----------
975	c Dust devil :
976	c ----------
977	IF(callddevil) then
978	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
979	& zdqdev,zdqsdev)
980
981	if (dustbin.ge.1) then
982	do iq=1,nq
983	DO l=1,nlayer
984	DO ig=1,ngrid
985	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
986	ENDDO
987	ENDDO
988	enddo
989	do iq=1,nq
990	DO ig=1,ngrid
991	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
992	ENDDO
993	enddo
994	endif ! of if (dustbin.ge.1)
995
996	END IF ! of IF (callddevil)
997
998	c -------------
999	c Sedimentation : acts also on water ice
1000	c -------------
1001	IF (sedimentation) THEN
1002	!call zerophys(ngridnlayernq, zdqsed)
1003	zdqsed(1:ngrid,1:nlayer,1:nq)=0
1004	!call zerophys(ngrid*nq, zdqssed)
1005	zdqssed(1:ngrid,1:nq)=0
1006
1007	call callsedim(ngrid,nlayer, ptimestep,
1008	& pplev,zzlev, pt, rdust, rice,
1009	& pq, pdq, zdqsed, zdqssed,nq)
1010
1011	DO iq=1, nq
1012	DO l=1,nlayer
1013	DO ig=1,ngrid
1014	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1015	ENDDO
1016	ENDDO
1017	ENDDO
1018	DO iq=1, nq
1019	DO ig=1,ngrid
1020	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1021	ENDDO
1022	ENDDO
1023	END IF ! of IF (sedimentation)
1024
1025	c 7d. Updates
1026	c ---------
1027
1028	DO iq=1, nq
1029	DO ig=1,ngrid
1030
1031	c ---------------------------------
1032	c Updating tracer budget on surface
1033	c ---------------------------------
1034	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1035
1036	ENDDO ! (ig)
1037	ENDDO ! (iq)
1038
1039	endif ! of if (tracer)
1040
1041	#ifndef MESOSCALE
1042	c-----------------------------------------------------------------------
1043	c 8. THERMOSPHERE CALCULATION
1044	c-----------------------------------------------------------------------
1045
1046	if (callthermos) then
1047	call thermosphere(pplev,pplay,dist_sol,
1048	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1049	& pt,pq,pu,pv,pdt,pdq,
1050	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1051
1052	DO l=1,nlayer
1053	DO ig=1,ngrid
1054	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1055	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1056	& +zdteuv(ig,l)
1057	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1058	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1059	DO iq=1, nq
1060	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1061	ENDDO
1062	ENDDO
1063	ENDDO
1064
1065	endif ! of if (callthermos)
1066	#endif
1067
1068	c-----------------------------------------------------------------------
1069	c 9. Surface and sub-surface soil temperature
1070	c-----------------------------------------------------------------------
1071	c
1072	c
1073	c 9.1 Increment Surface temperature:
1074	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1075
1076	DO ig=1,ngrid
1077	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1078	ENDDO
1079
1080	c Prescribe a cold trap at south pole (except at high obliquity !!)
1081	c Temperature at the surface is set there to be the temperature
1082	c corresponding to equilibrium temperature between phases of CO2
1083
1084	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
1085	#ifndef MESOSCALE
1086	if (caps.and.(obliquit.lt.27.)) then
1087	! NB: Updated surface pressure, at grid point 'ngrid', is
1088	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1089	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1090	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1091	endif
1092	#endif
1093	c -------------------------------------------------------------
1094	c Change of surface albedo (set to 0.4) in case of ground frost
1095	c everywhere except on the north permanent cap and in regions
1096	c covered by dry ice.
1097	c ALWAYS PLACE these lines after newcondens !!!
1098	c -------------------------------------------------------------
1099	do ig=1,ngrid
1100	if ((co2ice(ig).eq.0).and.
1101	& (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
1102	albedo(ig,1) = alb_surfice
1103	albedo(ig,2) = alb_surfice
1104	endif
1105	enddo ! of do ig=1,ngrid
1106	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
1107
1108	c
1109	c 9.2 Compute soil temperatures and subsurface heat flux:
1110	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1111	IF (callsoil) THEN
1112	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1113	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
1114	ENDIF
1115
1116	c-----------------------------------------------------------------------
1117	c 10. Write output files
1118	c ----------------------
1119
1120	c -------------------------------
1121	c Dynamical fields incrementation
1122	c -------------------------------
1123	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1124	! temperature, zonal and meridional wind
1125	DO l=1,nlayer
1126	DO ig=1,ngrid
1127	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1128	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1129	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1130	ENDDO
1131	ENDDO
1132
1133	! tracers
1134	DO iq=1, nq
1135	DO l=1,nlayer
1136	DO ig=1,ngrid
1137	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1138	ENDDO
1139	ENDDO
1140	ENDDO
1141
1142	! surface pressure
1143	DO ig=1,ngrid
1144	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1145	ENDDO
1146
1147	! pressure
1148	DO l=1,nlayer
1149	DO ig=1,ngrid
1150	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1151	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1152	ENDDO
1153	ENDDO
1154
1155	! Density
1156	DO l=1,nlayer
1157	DO ig=1,ngrid
1158	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1159	ENDDO
1160	ENDDO
1161
1162	! Potential Temperature
1163
1164	DO ig=1,ngridmx
1165	DO l=1,nlayermx
1166	zh(ig,l) = zt(ig,l)(zplay(ig,l)/zplev(ig,1))*rcp
1167	ENDDO
1168	ENDDO
1169
1170
1171	c Compute surface stress : (NB: z0 is a common in surfdat.h)
1172	c DO ig=1,ngrid
1173	c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
1174	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1175	c ENDDO
1176
1177	c Sum of fluxes in solar spectral bands (for output only)
1178	DO ig=1,ngrid
1179	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1180	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1181	ENDDO
1182	c ***** TEST ****************************************************
1183	ztim1 = 999
1184	DO l=1,nlayer
1185	DO ig=1,ngrid
1186	if (pt(ig,l).lt.ztim1) then
1187	ztim1 = pt(ig,l)
1188	igmin = ig
1189	lmin = l
1190	end if
1191	ENDDO
1192	ENDDO
1193	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1194	write(,) 'PHYSIQ: stability WARNING :'
1195	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1196	& 'ig l =', igmin, lmin
1197	end if
1198	c *******************************************************************
1199
1200	c ---------------------
1201	c Outputs to the screen
1202	c ---------------------
1203
1204	IF (lwrite) THEN
1205	PRINT*,'Global diagnostics for the physics'
1206	PRINT,'Variables and their increments x and dx/dt dt'
1207	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1208	WRITE(*,'(2f10.5,2f15.5)')
1209	s tsurf(igout),zdtsurf(igout)*ptimestep,
1210	s pplev(igout,1),pdpsrf(igout)*ptimestep
1211	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1212	WRITE(*,'(i4,6f10.5)') (l,
1213	s pu(igout,l),pdu(igout,l)*ptimestep,
1214	s pv(igout,l),pdv(igout,l)*ptimestep,
1215	s pt(igout,l),pdt(igout,l)*ptimestep,
1216	s l=1,nlayer)
1217	ENDIF ! of IF (lwrite)
1218
1219	IF (ngrid.NE.1) THEN
1220
1221	#ifndef MESOSCALE
1222	c -------------------------------------------------------------------
1223	c Writing NetCDF file "RESTARTFI" at the end of the run
1224	c -------------------------------------------------------------------
1225	c Note: 'restartfi' is stored just before dynamics are stored
1226	c in 'restart'. Between now and the writting of 'restart',
1227	c there will have been the itau=itau+1 instruction and
1228	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1229	c thus we store for time=time+dtvr
1230
1231	IF(lastcall) THEN
1232	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1233	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1234	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1235	. ptimestep,pday,
1236	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
1237	. area,albedodat,inertiedat,zmea,zstd,zsig,
1238	. zgam,zthe)
1239	ENDIF
1240	#endif
1241
1242	c -------------------------------------------------------------------
1243	c Calculation of diagnostic variables written in both stats and
1244	c diagfi files
1245	c -------------------------------------------------------------------
1246
1247	if (tracer) then
1248	if (water) then
1249
1250	call zerophys(ngrid,mtot)
1251	call zerophys(ngrid,icetot)
1252	call zerophys(ngrid,rave)
1253	call zerophys(ngrid,tauTES)
1254	do ig=1,ngrid
1255	do l=1,nlayermx
1256	mtot(ig) = mtot(ig) +
1257	& zq(ig,l,igcm_h2o_vap) *
1258	& (pplev(ig,l) - pplev(ig,l+1)) / g
1259	icetot(ig) = icetot(ig) +
1260	& zq(ig,l,igcm_h2o_ice) *
1261	& (pplev(ig,l) - pplev(ig,l+1)) / g
1262	rave(ig) = rave(ig) +
1263	& zq(ig,l,igcm_h2o_ice) *
1264	& (pplev(ig,l) - pplev(ig,l+1)) / g *
1265	& rice(ig,l) * (1.+nuice_ref)
1266	c Computing abs optical depth at 825 cm-1 in each
1267	c layer to simulate NEW TES retrieval
1268	Qabsice = min(
1269	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
1270	& )
1271	opTES(ig,l)= 0.75 * Qabsice *
1272	& zq(ig,l,igcm_h2o_ice) *
1273	& (pplev(ig,l) - pplev(ig,l+1)) / g
1274	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1275	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1276	enddo
1277	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
1278	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1279	enddo
1280
1281	endif ! of if (water)
1282	endif ! of if (tracer)
1283
1284	c -----------------------------------------------------------------
1285	c WSTATS: Saving statistics
1286	c -----------------------------------------------------------------
1287	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1288	c which can later be used to make the statistic files of the run:
1289	c "stats") only possible in 3D runs !
1290
1291	IF (callstats) THEN
1292
1293	call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1294	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1295	call wstats(ngrid,"co2ice","CO2 ice cover",
1296	& "kg.m-2",2,co2ice)
1297	call wstats(ngrid,"fluxsurf_lw",
1298	& "Thermal IR radiative flux to surface","W.m-2",2,
1299	& fluxsurf_lw)
1300	call wstats(ngrid,"fluxsurf_sw",
1301	& "Solar radiative flux to surface","W.m-2",2,
1302	& fluxsurf_sw_tot)
1303	call wstats(ngrid,"fluxtop_lw",
1304	& "Thermal IR radiative flux to space","W.m-2",2,
1305	& fluxtop_lw)
1306	call wstats(ngrid,"fluxtop_sw",
1307	& "Solar radiative flux to space","W.m-2",2,
1308	& fluxtop_sw_tot)
1309	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1310	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1311	call wstats(ngrid,"v","Meridional (North-South) wind",
1312	& "m.s-1",3,zv)
1313	call wstats(ngrid,"w","Vertical (down-up) wind",
1314	& "m.s-1",3,pw)
1315	call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
1316	call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
1317	c call wstats(ngrid,"q2",
1318	c & "Boundary layer eddy kinetic energy",
1319	c & "m2.s-2",3,q2)
1320	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1321	c & emis)
1322	c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
1323	c & 2,zstress)
1324	c call wstats(ngrid,"sw_htrt","sw heat.rate",
1325	c & "W.m-2",3,zdtsw)
1326	c call wstats(ngrid,"lw_htrt","lw heat.rate",
1327	c & "W.m-2",3,zdtlw)
1328
1329	if (tracer) then
1330	if (water) then
1331	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1332	& *mugaz/mmol(igcm_h2o_vap)
1333	call wstats(ngrid,"vmr_h2ovapor",
1334	& "H2O vapor volume mixing ratio","mol/mol",
1335	& 3,vmr)
1336	vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1337	& *mugaz/mmol(igcm_h2o_ice)
1338	call wstats(ngrid,"vmr_h2oice",
1339	& "H2O ice volume mixing ratio","mol/mol",
1340	& 3,vmr)
1341	call wstats(ngrid,"h2o_ice_s",
1342	& "surface h2o_ice","kg/m2",
1343	& 2,qsurf(1,igcm_h2o_ice))
1344
1345	call wstats(ngrid,"mtot",
1346	& "total mass of water vapor","kg/m2",
1347	& 2,mtot)
1348	call wstats(ngrid,"icetot",
1349	& "total mass of water ice","kg/m2",
1350	& 2,icetot)
1351	call wstats(ngrid,"reffice",
1352	& "Mean reff","m",
1353	& 2,rave)
1354	c call wstats(ngrid,"rice",
1355	c & "Ice particle size","m",
1356	c & 3,rice)
1357	c If activice is true, tauTES is computed in aeropacity.F;
1358	if (.not.activice) then
1359	call wstats(ngrid,"tauTESap",
1360	& "tau abs 825 cm-1","",
1361	& 2,tauTES)
1362	endif
1363
1364	endif ! of if (water)
1365
1366	if (thermochem.or.photochem) then
1367	do iq=1,nq
1368	if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
1369	. (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
1370	. (noms(iq).eq."h2").or.
1371	. (noms(iq).eq."o3")) then
1372	do l=1,nlayer
1373	do ig=1,ngrid
1374	vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
1375	end do
1376	end do
1377	call wstats(ngrid,"vmr_"//trim(noms(iq)),
1378	. "Volume mixing ratio","mol/mol",3,vmr)
1379	endif
1380	enddo
1381	endif ! of if (thermochem.or.photochem)
1382
1383	endif ! of if (tracer)
1384
1385	IF(lastcall) THEN
1386	write (,) "Writing stats..."
1387	call mkstats(ierr)
1388	ENDIF
1389
1390	ENDIF !if callstats
1391
1392	c (Store EOF for Mars Climate database software)
1393	IF (calleofdump) THEN
1394	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1395	ENDIF
1396
1397
1398	#ifdef MESOSCALE
1399	!!!
1400	!!! OUTPUT FIELDS
1401	!!!
1402	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
1403	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
1404	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
1405	IF (igcm_h2o_ice .ne. 0) THEN
1406	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
1407	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1408	. * mugaz / mmol(igcm_h2o_ice)
1409	ENDIF
1410	c AUTOMATICALLY GENERATED FROM REGISTRY
1411	#include "fill_save.inc"
1412	#else
1413
1414	c ==========================================================
1415	c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
1416	c any variable for diagnostic (output with period
1417	c "ecritphy", set in "run.def")
1418	c ==========================================================
1419	c WRITEDIAGFI can ALSO be called from any other subroutines
1420	c for any variables !!
1421	c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1422	c & emis)
1423	! call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
1424	! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
1425	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1426	& tsurf)
1427	call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1428	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
1429	& co2ice)
1430	c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
1431	c & "K",2,zt(1,7))
1432	call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1433	& fluxsurf_lw)
1434	call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1435	& fluxsurf_sw_tot)
1436	call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1437	& fluxtop_lw)
1438	call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1439	& fluxtop_sw_tot)
1440	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1441	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1442	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1443	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1444	call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1445	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1446	! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1447	c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
1448	c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
1449	c & zstress)
1450	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
1451	c & 'w.m-2',3,zdtsw)
1452	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
1453	c & 'w.m-2',3,zdtlw)
1454	c CALL WRITEDIAGFI(ngridmx,'tauTESap',
1455	c & 'tau abs 825 cm-1',
1456	c & '',2,tauTES)
1457
1458
1459	c ----------------------------------------------------------
1460	c Outputs of the CO2 cycle
1461	c ----------------------------------------------------------
1462
1463	if (tracer.and.(igcm_co2.ne.0)) then
1464	! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
1465	! & "kg/kg",2,zq(1,1,igcm_co2))
1466	! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
1467	! & "kg/kg",3,zq(1,1,igcm_co2))
1468
1469	! Compute co2 column
1470	call zerophys(ngrid,co2col)
1471	do l=1,nlayermx
1472	do ig=1,ngrid
1473	co2col(ig)=co2col(ig)+
1474	& zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
1475	enddo
1476	enddo
1477	call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
1478	& co2col)
1479	endif ! of if (tracer.and.(igcm_co2.ne.0))
1480
1481	c ----------------------------------------------------------
1482	c Outputs of the water cycle
1483	c ----------------------------------------------------------
1484	if (tracer) then
1485	if (water) then
1486
1487	CALL WRITEDIAGFI(ngridmx,'mtot',
1488	& 'total mass of water vapor',
1489	& 'kg/m2',2,mtot)
1490	CALL WRITEDIAGFI(ngridmx,'icetot',
1491	& 'total mass of water ice',
1492	& 'kg/m2',2,icetot)
1493	c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1494	c & *mugaz/mmol(igcm_h2o_ice)
1495	c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
1496	c & 'mol/mol',3,vmr)
1497	CALL WRITEDIAGFI(ngridmx,'reffice',
1498	& 'Mean reff',
1499	& 'm',2,rave)
1500	c call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
1501	c & 'm',3,rice)
1502	c If activice is true, tauTES is computed in aeropacity.F;
1503	if (.not.activice) then
1504	CALL WRITEDIAGFI(ngridmx,'tauTESap',
1505	& 'tau abs 825 cm-1',
1506	& '',2,tauTES)
1507	endif
1508	call WRITEDIAGFI(ngridmx,'h2o_ice_s',
1509	& 'surface h2o_ice',
1510	& 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
1511	endif !(water)
1512
1513
1514	if (water.and..not.photochem) then
1515	iq=nq
1516	c write(str2(1:2),'(i2.2)') iq
1517	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
1518	c & 'kg.m-2',2,zdqscloud(1,iq))
1519	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
1520	c & 'kg/kg',3,zdqchim(1,1,iq))
1521	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
1522	c & 'kg/kg',3,zdqdif(1,1,iq))
1523	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
1524	c & 'kg/kg',3,zdqadj(1,1,iq))
1525	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
1526	c & 'kg/kg',3,zdqc(1,1,iq))
1527	endif !(water.and..not.photochem)
1528	endif
1529
1530	c ----------------------------------------------------------
1531	c Outputs of the dust cycle
1532	c ----------------------------------------------------------
1533
1534	c call WRITEDIAGFI(ngridmx,'tauref',
1535	c & 'Dust ref opt depth','NU',2,tauref)
1536
1537	if (tracer.and.(dustbin.ne.0)) then
1538	c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
1539	if (doubleq) then
1540	c call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
1541	c & 'kg.m-2',2,qsurf(1,1))
1542	c call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
1543	c & 'N.m-2',2,qsurf(1,2))
1544	c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
1545	c & 'kg.m-2.s-1',2,zdqsdev(1,1))
1546	c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
1547	c & 'kg.m-2.s-1',2,zdqssed(1,1))
1548	call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
1549	& 'm',3,rdust*ref_r0)
1550	call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
1551	& 'kg/kg',3,pq(1,1,igcm_dust_mass))
1552	c call WRITEDIAGFI(ngridmx,'dustN','Dust number',
1553	c & 'part/kg',3,pq(1,1,igcm_dust_number))
1554	else
1555	do iq=1,dustbin
1556	write(str2(1:2),'(i2.2)') iq
1557	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
1558	& 'kg/kg',3,zq(1,1,iq))
1559	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
1560	& 'kg.m-2',2,qsurf(1,iq))
1561	end do
1562	endif ! (doubleq)
1563	c if (submicron) then
1564	c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
1565	c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))
1566	c endif ! (submicron)
1567	end if ! (tracer.and.(dustbin.ne.0))
1568
1569	c ----------------------------------------------------------
1570	c ----------------------------------------------------------
1571	c PBL OUTPUS
1572	c ----------------------------------------------------------
1573	c ----------------------------------------------------------
1574
1575
1576	c ----------------------------------------------------------
1577	c Outputs of surface layer
1578	c ----------------------------------------------------------
1579
1580
1581	z_out=0.
1582	if (calltherm .and. (z_out .gt. 0.)) then
1583	call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
1584	& ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar)
1585
1586	zu2(:)=sqrt(zu(:,1)zu(:,1)+zv(:,1)zv(:,1))
1587	call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
1588	& 'horizontal velocity norm','m/s',
1589	& 2,zu2)
1590
1591	call WRITEDIAGFI(ngridmx,'Teta_out',
1592	& 'potential temperature at z_out','K',
1593	& 2,Teta_out)
1594	call WRITEDIAGFI(ngridmx,'u_out',
1595	& 'horizontal velocity norm at z_out','m/s',
1596	& 2,u_out)
1597	call WRITEDIAGFI(ngridmx,'u*',
1598	& 'friction velocity','m/s',
1599	& 2,ustar)
1600	call WRITEDIAGFI(ngridmx,'teta*',
1601	& 'friction potential temperature','K',
1602	& 2,tstar)
1603	else
1604	if((.not. calltherm).and.(z_out .gt. 0.)) then
1605	print*, 'WARNING : no interpolation in surface-layer :'
1606	print*, 'Outputing surface-layer quantities without thermals
1607	& does not make sense'
1608	endif
1609	endif
1610
1611	c ----------------------------------------------------------
1612	c Outputs of thermals
1613	c ----------------------------------------------------------
1614	if (calltherm) then
1615
1616	! call WRITEDIAGFI(ngrid,'dtke',
1617	! & 'tendance tke thermiques','m2/s2',
1618	! & 3,dtke_th)
1619	! call WRITEDIAGFI(ngrid,'d_u_ajs',
1620	! & 'tendance u thermiques','m/s',
1621	! & 3,pdu_th*ptimestep)
1622	! call WRITEDIAGFI(ngrid,'d_v_ajs',
1623	! & 'tendance v thermiques','m/s',
1624	! & 3,pdv_th*ptimestep)
1625	! if (tracer) then
1626	! if (nq .eq. 2) then
1627	! call WRITEDIAGFI(ngrid,'deltaq_th',
1628	! & 'delta q thermiques','kg/kg',
1629	! & 3,ptimestep*pdq_th(:,:,2))
1630	! endif
1631	! endif
1632
1633	lmax_th_out(:)=real(lmax_th(:))
1634
1635	call WRITEDIAGFI(ngridmx,'lmax_th',
1636	& 'hauteur du thermique','K',
1637	& 2,lmax_th_out)
1638	call WRITEDIAGFI(ngridmx,'hfmax_th',
1639	& 'maximum TH heat flux','K.m/s',
1640	& 2,hfmax_th)
1641	call WRITEDIAGFI(ngridmx,'wmax_th',
1642	& 'maximum TH vertical velocity','m/s',
1643	& 2,wmax_th)
1644
1645	endif
1646
1647	c ----------------------------------------------------------
1648	c ----------------------------------------------------------
1649	c END OF PBL OUTPUS
1650	c ----------------------------------------------------------
1651	c ----------------------------------------------------------
1652
1653
1654	c ----------------------------------------------------------
1655	c Output in netcdf file "diagsoil.nc" for subterranean
1656	c variables (output every "ecritphy", as for writediagfi)
1657	c ----------------------------------------------------------
1658
1659	! Write soil temperature
1660	! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
1661	! & 3,tsoil)
1662	! Write surface temperature
1663	! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
1664	! & 2,tsurf)
1665
1666	c ==========================================================
1667	c END OF WRITEDIAGFI
1668	c ==========================================================
1669	#endif
1670
1671	ELSE ! if(ngrid.eq.1)
1672
1673	print,'Ls =',zls180./pi,
1674	& ' tauref(700 Pa) =',tauref
1675	c ----------------------------------------------------------------------
1676	c Output in grads file "g1d" (ONLY when using testphys1d)
1677	c (output at every X physical timestep)
1678	c ----------------------------------------------------------------------
1679	c
1680	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
1681	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
1682	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
1683
1684	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
1685	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
1686	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
1687	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
1688
1689	! THERMALS STUFF 1D
1690
1691	z_out=0.
1692	if (calltherm .and. (z_out .gt. 0.)) then
1693	call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
1694	& ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar)
1695
1696	zu2(:)=sqrt(zu(:,1)zu(:,1)+zv(:,1)zv(:,1))
1697	call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
1698	& 'horizontal velocity norm','m/s',
1699	& 0,zu2)
1700
1701	call WRITEDIAGFI(ngridmx,'Teta_out',
1702	& 'potential temperature at z_out','K',
1703	& 0,Teta_out)
1704	call WRITEDIAGFI(ngridmx,'u_out',
1705	& 'horizontal velocity norm at z_out','m/s',
1706	& 0,u_out)
1707	call WRITEDIAGFI(ngridmx,'u*',
1708	& 'friction velocity','m/s',
1709	& 0,ustar)
1710	call WRITEDIAGFI(ngridmx,'teta*',
1711	& 'friction potential temperature','K',
1712	& 0,tstar)
1713	else
1714	if((.not. calltherm).and.(z_out .gt. 0.)) then
1715	print*, 'WARNING : no interpolation in surface-layer :'
1716	print*, 'Outputing surface-layer quantities without thermals
1717	& does not make sense'
1718	endif
1719	endif
1720
1721	! ---
1722
1723
1724
1725	if(calltherm) then
1726
1727	lmax_th_out(:)=real(lmax_th(:))
1728
1729	call WRITEDIAGFI(ngridmx,'lmax_th',
1730	& 'hauteur du thermique','point',
1731	& 0,lmax_th_out)
1732	call WRITEDIAGFI(ngridmx,'hfmax_th',
1733	& 'maximum TH heat flux','K.m/s',
1734	& 0,hfmax_th)
1735	call WRITEDIAGFI(ngridmx,'wmax_th',
1736	& 'maximum TH vertical velocity','m/s',
1737	& 0,wmax_th)
1738
1739	co2col(:)=0.
1740	if (tracer) then
1741	do l=1,nlayermx
1742	do ig=1,ngrid
1743	co2col(ig)=co2col(ig)+
1744	& zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
1745	enddo
1746	enddo
1747
1748	end if
1749	call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' &
1750	& ,'kg/m-2',0,co2col)
1751	endif
1752	call WRITEDIAGFI(ngrid,'w','vertical velocity' &
1753	& ,'m/s',1,pw)
1754	call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
1755	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
1756	& tsurf)
1757	! call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
1758	! call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
1759
1760	call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
1761	call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
1762	! or output in diagfi.nc (for testphys1d)
1763	call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
1764	call WRITEDIAGFI(ngridmx,'temp','Temperature',
1765	& 'K',1,zt)
1766
1767	if(tracer) then
1768	c CALL writeg1d(ngrid,1,tau,'tau','SI')
1769	do iq=1,nq
1770	c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
1771	call WRITEDIAGFI(ngridmx,trim(noms(iq)),
1772	& trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
1773	end do
1774	end if
1775
1776	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
1777
1778	do l=2,nlayer-1
1779	tmean=zt(1,l)
1780	if(zt(1,l).ne.zt(1,l-1))
1781	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
1782	zlocal(l)= zlocal(l-1)
1783	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
1784	enddo
1785	zlocal(nlayer)= zlocal(nlayer-1)-
1786	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
1787	& rnew(1,nlayer)*tmean/g
1788
1789	END IF ! if(ngrid.ne.1)
1790
1791	icount=icount+1
1792	RETURN
1793	END

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

Download in other formats:

Original Format