Context Navigation

physiq.F @ 162

Last change on this file since 162 was 161, checked in by acolaitis, 14 years ago

================================================
======== IMPLEMENTATION OF THERMALS ============
================================================

Author: A. Colaitis (2011-06-16)

The main goal of this revision is to start including the thermals into the model
for development purposes. Users should not use the thermals yet, as
several major configuration changes still need to be done.

This version includes :

updraft and downdraft parametrizations
velocity in the thermal, including drag
plume height analysis
closure equation
updraft transport of heat, tracers and momentum
downdraft transport of heat

This model should not be used without upcoming developments, namely :

downdraft transport of tracers and momentum
updraft & downdraft transport of q2 (tke)
revision of vdif_kc to compute q2 for non-stratified cases

Thermals could also include in a later revision :

momentum loss during transport (horizontal drag)

Compilation of the thermals has been successfully tested on ifort, gfortran and pgf90

================================================
================================================

M libf/phymars/callkeys.h
M libf/phymars/inifis.F

Added new control flags to call the thermals :

calltherm (false by default) <- to call thermals
outptherm (false by default) <- to output thermal-related diagnostics (for dev purposes)

================================================

M libf/phymars/vdifc.F
^{------> added a temporary output for thermal-related diagnostics}

M libf/phymars/testphys1d.F
^{------> added treatment for a initialization from a profile of neutral gas (ar)}

-> will be transformed in a decaying tracer for thermal diagnostics

M libf/phymars/physiq.F
^{------> added a section to call the thermals}

-> changed the call to convadj
-> added thermal-related outputs for diagnostics

M libf/phymars/convadj.F
^{------> takes now into account the height of thermals to execute convective adjustment}

=> note : convective adjustment needs to be activated when using thermals, in case of a

second instable layer above the thermals

================================================

A libf/phymars/calltherm_interface.F90
^{------> Interface between physiq.F and the thermals}

A libf/phymars/calltherm_mars.F90
^{------> Routine running the sub-timestep of the thermals}

A libf/phymars/thermcell_main_mars.F90
^{------> Main thermals routine specific to Martian physics}

A libf/phymars/thermcell_dqupdown.F90
^{------> Thermals subroutine computing transport of quantities by updrafts and downdrafts}

A libf/phymars/thermcell.F90
^{------> Module including parameters from the Earth to Mars importation. Will disappear in future dev}

================================================
================================================

File size: 56.8 KB

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

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

Download in other formats:

Original Format