Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

physiq.F @ 54

Last change on this file since 54 was 54, checked in by aslmd, 14 years ago

LMD_LES_MARS: Tests with new physics in unified frame
M 53 mesoscale/LMD_LES_MARS/modif_mars/module_model_constants.F
M 53 mesoscale/LMD_LES_MARS/modif_mars/module_first_rk_step_part1.F
M 53 mesoscale/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F
Modifications made to cope with new soil scheme

M 53 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_physiq.F
Commented aeronomars in old physics so that the model can compile with lots of grid points

LMD_MM_MARS: Implement tracers with the new physics. Test with JBM runs with 5 tracers

including radiatively active dust and water ice

TODO: propagate reading "traceur.def" and few initialisations from

testphys1d.F into module_lmd_driver.F

M 53 mars/libf/phymars/dimradmars.h
M 53 mars/libf/phymars/callradite.F
Change it back to both dust and ice radiatively active

M 53 mesoscale/LMD_MM_MARS/makemeso
Allow to "fresh start" with the inclusion of "clean -a" in the script

M 53 mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM
M 53 mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F
M 53 mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/solve_em.F
M 53 mesoscale/LMD_MM_MARS/SRC/WRFV2/main/real_em.F
M 53 mesoscale/LMD_MM_MARS/SIMU/runmeso
Added the option mars==11 with 5 tracers: CO2 H2Ovap H2Oice DUST DUSTN
-- Note that the order matters for H2Ovap and H2Oice

LMD_MM_MARS: Create a folder to run LMD GCM with exact same physics as mesoscale

A 0 mesoscale/LMDZ.MARS/libf/phymars/physiq.F
All files in mesoscale/LMDZ.MARS/libf are links except this one
--> Because specific WRITEDIAGFI commands are set to output what is needed PREP_MARS
--> Something must be done so that qCO2, qdust, qdustN can be propagated too !!!

A 0 mesoscale/LMDZ.MARS/libf_gcm
A 0 mesoscale/LMDZ.MARS/in_lmdz_mars_newphys
Two important links that must not be broken

A 0 mesoscale/LMDZ.MARS/myGCM
This folder is for runs. Only links here -- except temp files which are not synchronized.

M 53 mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/compile
A 0 mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/run.def
A 0 mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/traceur.def
M 53 mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/callphys.def
A 0 mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/callphys.def.csttau
Those files are needed to run the GCM in order to prepare initial & boundary

conditions for the mesoscale... but basically it is the exact same GCM

--> compile is a convenient script

A 0 mesoscale/LMDZ.MARS/makegcm
A 0 mesoscale/LMDZ.MARS/makegcm_g95
A 0 mesoscale/LMDZ.MARS/create_make_gcm
This is supposed to be merged one day with files in trunk/mars/

A 0 mesoscale/LMDZ.MARS/myGCM/DEFS_JB
Necessary files to restart a GCM run corresponding to new water cycle in JB's thesis

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

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

Download in other formats:

Original Format