Context Navigation

physiq.F @ 411

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

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

Download in other formats:

Original Format