Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

physiq.F_old @ 1524

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

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

Download in other formats:

Original Format