Context Navigation

physiq.F @ 190

Last change on this file since 190 was 185, checked in by acolaitis, 14 years ago

17/06/2011 == AC

Added new settings for the Martian thermals from new LES observations
Revamped thermcell's module variables to allow it's removal
Minor changes in physiq and meso_physiq for the call to thermals
Switched from dynamic to static memory allocation for all thermals variable

to gain computation speed

Thermals now output maximum speed, maximum heat flux, and maximum height in thermal plumes

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

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

Download in other formats:

Original Format