Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

physiq.F @ 1647

Last change on this file since 1647 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 46.4 KB

Line
1	SUBROUTINE physiq(ngrid,nlayer,nq,
2	$ firstcall,lastcall,
3	$ pday,ptime,ptimestep,
4	$ pplev,pplay,pphi,
5	$ pu,pv,pt,pq,
6	$ pw,
7	$ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
8
9
10	IMPLICIT NONE
11	c=======================================================================
12	c
13	c subject:
14	c --------
15	c
16	c Organisation of the physical parametrisations of the LMD
17	c martian atmospheric general circulation model.
18	c
19	c The GCM can be run without or with tracer transport
20	c depending on the value of Logical "tracer" in file "callphys.def"
21	c Tracers may be water vapor, ice OR chemical species OR dust particles
22	c
23	c SEE comments in initracer.F about numbering of tracer species...
24	c
25	c It includes:
26	c
27	c 1. Initialization:
28	c 1.1 First call initializations
29	c 1.2 Initialization for every call to physiq
30	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
31	c 2. Compute radiative transfer tendencies
32	c (longwave and shortwave) for CO2 and dust.
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 10. 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
61	c
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) /
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 "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 "fisice.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	SAVE day_ini, icount
185	SAVE aerosol, tsurf,tsoil
186	SAVE co2ice,albedo,emis, q2
187	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
188	SAVE ig_vl1
189
190	REAL stephan
191	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
192	SAVE stephan
193
194	c Local variables :
195	c -----------------
196
197
198	CHARACTER*80 fichier
199	INTEGER l,ig,ierr,igout,iq, tapphys
200	INTEGER iqmin ! Used if iceparty engaged
201
202	REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)
203	REAL fluxsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)
204	REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)
205	REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)
206	c for clear area (uncovered by clouds) :
207	REAL clsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)
208	REAL clsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)
209	REAL cltop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)
210	REAL cltop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)
211	REAL tauref(ngridmx) ! Reference column optical depth at 700 Pa
212	! (used if active=F)
213	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
214	REAL zls ! solar longitude (rad)
215	REAL zday ! date (time since Ls=0, in martian days)
216	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
217	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
218	REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
219
220
221	c Tendancies due to various processes:
222	REAL dqsurf(ngridmx,nqmx)
223	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
224	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
225	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear area
226	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear area
227	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
228	REAL zdtnlte(ngridmx,nlayermx) ! (K/s)
229	REAL zdtsurf(ngridmx) ! (K/s)
230	REAL zdtcloud(ngridmx,nlayermx)
231	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
232	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
233	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
234	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
235	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
236	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
237	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
238	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
239
240	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
241	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
242	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
243	REAL zdqadj(ngridmx,nlayermx,nqmx)
244	REAL zdqc(ngridmx,nlayermx,nqmx)
245	REAL zdqcloud(ngridmx,nlayermx,nqmx)
246	REAL zdqscloud(ngridmx,nqmx)
247	REAL zdqchim(ngridmx,nlayermx,nqmx)
248	REAL zdqschim(ngridmx,nqmx)
249
250	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
251	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
252	REAL zdumolvis(ngridmx,nlayermx)
253	REAL zdvmolvis(ngridmx,nlayermx)
254	real zdqmoldiff(ngridmx,nlayermx,nqmx)
255
256	c Local variable for local intermediate calcul:
257	REAL zflubid(ngridmx)
258	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
259	REAL zdum1(ngridmx,nlayermx)
260	REAL zdum2(ngridmx,nlayermx)
261	REAL ztim1,ztim2,ztim3, z1,z2
262	REAL ztime_fin
263	REAL zdh(ngridmx,nlayermx)
264	REAL pclc_min ! minimum of the cloud fraction over the whole domain
265	INTEGER length
266	PARAMETER (length=100)
267
268	c local variables only used for diagnostic (output in file "diagfi" or "stats")
269	c -----------------------------------------------------------------------------
270	REAL ps(ngridmx), zt(ngridmx,nlayermx)
271	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
272	REAL zq(ngridmx,nlayermx,nqmx)
273	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
274	character*2 str2
275	character*5 str5
276	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
277	real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
278	real qtot1,qtot2 ! total aerosol mass
279	integer igmin, lmin
280	logical tdiag
281
282
283	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
284	REAL zstress(ngrid), cd
285	real hco2(nqmx),tmean, zlocal(nlayermx)
286	real rho(ngridmx,nlayermx) ! density
287	real vmr(ngridmx,nlayermx) ! volume mixing ratio
288
289	REAL time_phys
290
291	c=======================================================================
292
293	c 1. Initialisation:
294	c -----------------
295
296	c 1.1 Initialisation only at first call
297	c ---------------------------------------
298	IF (firstcall) THEN
299
300	c variables set to 0
301	c ~~~~~~~~~~~~~~~~~~
302	call zerophys(ngridnlayernaerkind,aerosol)
303	call zerophys(ngrid*nlayer,dtrad)
304	call zerophys(ngrid,fluxrad)
305
306	c read startfi
307	c ~~~~~~~~~~~~
308
309	! Read netcdf initial physical parameters.
310	CALL phyetat0 ("startfi.nc",0,0,
311	& nsoilmx,nq,
312	& day_ini,time_phys,
313	& tsurf,tsoil,emis,q2,qsurf,co2ice)
314
315	if (pday.ne.day_ini) then
316	write(,) "PHYSIQ: ERROR: bad synchronization between ",
317	& "physics and dynamics"
318	write(,) "dynamics day: ",pday
319	write(,) "physics day: ",day_ini
320	stop
321	endif
322
323	write (,) 'In physiq day_ini =', day_ini
324
325	c Initialize albedo and orbital calculation
326	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
327	CALL surfini(ngrid,co2ice,qsurf,albedo)
328	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
329
330	c initialize soil
331	c ~~~~~~~~~~~~~~~
332	IF (callsoil) THEN
333	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
334	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
335	ELSE
336	PRINT*,'WARNING! Thermal conduction in the soil turned off'
337	DO ig=1,ngrid
338	capcal(ig)=1.e5
339	fluxgrd(ig)=0.
340	ENDDO
341	ENDIF
342	icount=1
343
344
345	c initialize tracers
346	c ~~~~~~~~~~~~~~~~~~
347	tracerdyn=tracer
348	IF (tracer) THEN
349	CALL initracer(qsurf,co2ice)
350	ENDIF ! end tracer
351
352	c Determining gridpoint near Viking Lander 1 (used for diagnostic only)
353	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
354
355	if(ngrid.ne.1) then
356	latvl1= 22.27
357	lonvl1= -47.94
358	ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim +
359	& int(1.5+(lonvl1+180)*iim/360.)
360	write(,) 'Viking Lander 1 GCM point: lat,lon',
361	& lati(ig_vl1)180/pi, long(ig_vl1)180/pi
362	end if
363
364	c Initialize thermospheric parameters
365	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
366
367	if (callthermos) call param_read
368
369	c Initialize R and Cp as constant
370
371	if (.not.callthermos .and. .not.photochem) then
372	do l=1,nlayermx
373	do ig=1,ngridmx
374	rnew(ig,l)=r
375	cpnew(ig,l)=cpp
376	mmean(ig,l)=mugaz
377	enddo
378	enddo
379	endif
380
381	ENDIF ! (end of "if firstcall")
382
383	c ---------------------------------------------------
384	c 1.2 Initializations done at every physical timestep:
385	c ---------------------------------------------------
386	c
387	IF (ngrid.NE.ngridmx) THEN
388	PRINT*,'STOP in PHYSIQ'
389	PRINT*,'Probleme de dimensions :'
390	PRINT*,'ngrid = ',ngrid
391	PRINT*,'ngridmx = ',ngridmx
392	STOP
393	ENDIF
394
395	c Initialize various variables
396	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
397	call zerophys(ngrid*nlayer, pdv)
398	call zerophys(ngrid*nlayer, pdu)
399	call zerophys(ngrid*nlayer, pdt)
400	call zerophys(ngridnlayernq, pdq)
401	call zerophys(ngrid, pdpsrf)
402	call zerophys(ngrid, zflubid)
403	call zerophys(ngrid, zdtsurf)
404	call zerophys(ngrid*nq, dqsurf)
405	igout=ngrid/2+1
406
407
408	zday=pday+ptime ! compute time, in sols (and fraction thereof)
409
410	c Compute Solar Longitude (Ls) :
411	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
412	if (season) then
413	call solarlong(zday,zls)
414	else
415	call solarlong(float(day_ini),zls)
416	end if
417
418	c Compute geopotential at interlayers
419	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
420	c ponderation des altitudes au niveau des couches en dp/p
421
422	DO l=1,nlayer
423	DO ig=1,ngrid
424	zzlay(ig,l)=pphi(ig,l)/g
425	ENDDO
426	ENDDO
427	DO ig=1,ngrid
428	zzlev(ig,1)=0.
429	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
430	ENDDO
431	DO l=2,nlayer
432	DO ig=1,ngrid
433	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
434	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
435	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
436	ENDDO
437	ENDDO
438
439
440	! Potential temperature calculation not the same in physiq and dynamic
441
442	c Compute potential temperature
443	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
444	DO l=1,nlayer
445	DO ig=1,ngrid
446	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
447	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
448	ENDDO
449	ENDDO
450
451	c-----------------------------------------------------------------------
452	c 1.2.5 Compute mean mass, cp, and R
453	c --------------------------------
454
455	if(photochem.or.callthermos) then
456	call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
457	endif
458	c-----------------------------------------------------------------------
459	c 2. Compute radiative tendencies :
460	c------------------------------------
461
462
463	IF (callrad) THEN
464	IF( MOD(icount-1,iradia).EQ.0) THEN
465
466	c Local Solar zenith angle
467	c ~~~~~~~~~~~~~~~~~~~~~~~~
468	CALL orbite(zls,dist_sol,declin)
469
470	IF(diurnal) THEN
471	ztim1=SIN(declin)
472	ztim2=COS(declin)COS(2.pi*(zday-.5))
473	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
474
475	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
476	s ztim1,ztim2,ztim3, mu0,fract)
477
478	ELSE
479	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
480	ENDIF
481
482	c NLTE cooling from CO2 emission
483	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
484
485	IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
486
487	c Find number of layers for LTE radiation calculations
488	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
489	& CALL nlthermeq(ngrid,nlayer,pplev,pplay)
490
491
492	c Atmospheric dust opacity and aerosol distribution:
493	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
494	CALL dustopacity(ngrid,nlayer,nq,zday,pplay,pplev,zls,pq,
495	$ tauref,tau,aerosol)
496
497	c Call main radiative transfer scheme
498	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
499	c Transfer through dust and CO2, except NIR CO2 absorption
500
501	CALL callradite(icount,ngrid,nlayer,aerosol,albedo,
502	$ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
503	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw)
504
505	c CO2 near infrared absorption
506	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
507	call zerophys(ngrid*nlayer,zdtnirco2)
508	if (callnirco2) then
509	call nirco2abs (ngrid,nlayer,pplay,dist_sol,
510	. mu0,fract,declin, zdtnirco2)
511	endif
512
513	c Radiative flux from the sky absorbed by the surface (W.m-2)
514	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
515	DO ig=1,ngrid
516	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
517	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
518	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
519	ENDDO
520
521
522	c Net atmospheric radiative heating rate (K.s-1)
523	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
524	IF(callnlte) THEN
525	CALL blendrad(ngrid, nlayer, pplay,
526	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
527	ELSE
528	DO l=1,nlayer
529	DO ig=1,ngrid
530	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
531	& +zdtnirco2(ig,l)
532	ENDDO
533	ENDDO
534	ENDIF
535
536
537
538	ENDIF ! of if(mod(icount-1,iradia).eq.0)
539
540
541
542	c Transformation of the radiative tendencies:
543	c -------------------------------------------
544
545	c Net radiative surface flux (W.m-2)
546	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
547	c
548	DO ig=1,ngrid
549	zplanck(ig)=tsurf(ig)*tsurf(ig)
550	zplanck(ig)=emis(ig)*
551	$ stephanzplanck(ig)zplanck(ig)
552	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
553	ENDDO
554
555
556	DO l=1,nlayer
557	DO ig=1,ngrid
558	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
559	ENDDO
560	ENDDO
561
562	ENDIF ! of IF (callrad)
563
564	c-----------------------------------------------------------------------
565	c 3. Gravity wave and subgrid scale topography drag :
566	c -------------------------------------------------
567
568
569	IF(calllott)THEN
570
571	CALL calldrag_noro(ngrid,nlayer,ptimestep,
572	& pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
573
574	DO l=1,nlayer
575	DO ig=1,ngrid
576	pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
577	pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
578	pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
579	ENDDO
580	ENDDO
581	ENDIF
582
583	c-----------------------------------------------------------------------
584	c 4. Vertical diffusion (turbulent mixing):
585	c -----------------------------------------
586	c
587	IF (calldifv) THEN
588
589
590	DO ig=1,ngrid
591	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
592	ENDDO
593
594	CALL zerophys(ngrid*nlayer,zdum1)
595	CALL zerophys(ngrid*nlayer,zdum2)
596	do l=1,nlayer
597	do ig=1,ngrid
598	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
599	enddo
600	enddo
601
602	c Calling vdif (Martian version WITH CO2 condensation)
603	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
604	$ ptimestep,capcal,lwrite,
605	$ pplay,pplev,zzlay,zzlev,z0,
606	$ pu,pv,zh,pq,tsurf,emis,qsurf,
607	$ zdum1,zdum2,zdh,pdq,zflubid,
608	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
609	& zdqdif,zdqsdif)
610
611	DO l=1,nlayer
612	DO ig=1,ngrid
613	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
614	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
615	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
616
617	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
618
619	ENDDO
620	ENDDO
621
622	DO ig=1,ngrid
623	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
624	ENDDO
625
626	if (tracer) then
627	DO iq=1, nq
628	DO l=1,nlayer
629	DO ig=1,ngrid
630	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
631	ENDDO
632	ENDDO
633	ENDDO
634	DO iq=1, nq
635	DO ig=1,ngrid
636	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
637	ENDDO
638	ENDDO
639	end if ! of if (tracer)
640
641	ELSE
642	DO ig=1,ngrid
643	zdtsurf(ig)=zdtsurf(ig)+
644	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
645	ENDDO
646	ENDIF ! of IF (calldifv)
647
648
649	c-----------------------------------------------------------------------
650	c 5. Dry convective adjustment:
651	c -----------------------------
652
653	IF(calladj) THEN
654
655	DO l=1,nlayer
656	DO ig=1,ngrid
657	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
658	ENDDO
659	ENDDO
660	CALL zerophys(ngrid*nlayer,zduadj)
661	CALL zerophys(ngrid*nlayer,zdvadj)
662	CALL zerophys(ngrid*nlayer,zdhadj)
663
664	CALL convadj(ngrid,nlayer,nq,ptimestep,
665	$ pplay,pplev,zpopsk,
666	$ pu,pv,zh,pq,
667	$ pdu,pdv,zdh,pdq,
668	$ zduadj,zdvadj,zdhadj,
669	$ zdqadj)
670
671	DO l=1,nlayer
672	DO ig=1,ngrid
673	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
674	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
675	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
676
677	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
678	ENDDO
679	ENDDO
680
681	if(tracer) then
682	DO iq=1, nq
683	DO l=1,nlayer
684	DO ig=1,ngrid
685	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
686	ENDDO
687	ENDDO
688	ENDDO
689	end if
690	ENDIF ! of IF(calladj)
691
692	c-----------------------------------------------------------------------
693	c 6. Carbon dioxide condensation-sublimation:
694	c -------------------------------------------
695
696	IF (callcond) THEN
697	CALL newcondens(ngrid,nlayer,nq,ptimestep,
698	$ capcal,pplay,pplev,tsurf,pt,
699	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
700	$ co2ice,albedo,emis,
701	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
702	$ fluxsurf_sw)
703
704	DO l=1,nlayer
705	DO ig=1,ngrid
706	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
707	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
708	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
709	ENDDO
710	ENDDO
711	DO ig=1,ngrid
712	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
713	ENDDO
714
715	IF (tracer) THEN
716	DO iq=1, nq
717	DO l=1,nlayer
718	DO ig=1,ngrid
719	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
720	ENDDO
721	ENDDO
722	ENDDO
723	ENDIF ! of IF (tracer)
724
725	ENDIF ! of IF (callcond)
726
727	c-----------------------------------------------------------------------
728	c 7. Specific parameterizations for tracers
729	c: -----------------------------------------
730
731	if (tracer) then
732
733	c 7a. Water and ice
734	c ---------------
735
736	c ---------------------------------------
737	c Water ice condensation in the atmosphere
738	c ----------------------------------------
739	IF (water) THEN
740	call watercloud(ngrid,nlayer, ptimestep,
741	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
742	& pq,pdq,zdqcloud,qsurf,zdqscloud,zdtcloud,
743	& nq,naerkind,tau,icount,zls)
744
745	if (activice) then
746	c Temperature variation due to latent heat release
747	DO l=1,nlayer
748	DO ig=1,ngrid
749	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
750	ENDDO
751	ENDDO
752	endif
753
754	IF (iceparty) then
755	iqmin=nq-1
756	ELSE
757	iqmin=nq
758	ENDIF
759
760	DO iq=iqmin,nq
761	DO l=1,nlayer
762	DO ig=1,ngrid
763	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq)
764	ENDDO
765	ENDDO
766	DO ig=1,ngrid
767	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq)
768	ENDDO
769	ENDDO
770
771	END IF ! of IF (water)
772
773	c 7b. Chemical species
774	c ------------------
775
776	c --------------
777	c photochemistry :
778	c --------------
779	IF (photochem .or. thermochem) then
780	call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
781	. zzlay,zday,pq,pdq,zdqchim,zdqschim,zdqcloud,zdqscloud)
782
783	c Photochemistry includes condensation of H2O2
784
785	do iq=nqchem_min,nq
786	if (noms(iq).eq."h2o2") then
787	DO l=1,nlayer
788	DO ig=1,ngrid
789	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq)
790	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq)
791	ENDDO
792	ENDDO
793	else
794	DO l=1,nlayer
795	DO ig=1,ngrid
796	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq)
797	ENDDO
798	ENDDO
799	endif
800	ENDDO
801	do iq=nqchem_min,nq
802	if (noms(iq).eq."h2o2") then
803	DO ig=1,ngrid
804	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq)
805	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq)
806	ENDDO
807	else
808	DO ig=1,ngrid
809	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq)
810	ENDDO
811	endif
812	ENDDO
813
814	END IF ! of IF (photochem.or.thermochem)
815
816	c 7c. Aerosol particles
817	c -------------------
818
819	c ----------
820	c Dust devil :
821	c ----------
822	IF(callddevil) then
823	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
824	& zdqdev,zdqsdev)
825
826	if (dustbin.ge.1) then
827	do iq=1,nq
828	DO l=1,nlayer
829	DO ig=1,ngrid
830	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
831	ENDDO
832	ENDDO
833	enddo
834	do iq=1,nq
835	DO ig=1,ngrid
836	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
837	ENDDO
838	enddo
839	endif ! of if (dustbin.ge.1)
840
841	END IF ! of IF (callddevil)
842
843	c -------------
844	c Sedimentation : acts also on water ice
845	c -------------
846	IF (sedimentation) THEN
847	call zerophys(ngridnlayernq, zdqsed)
848	call zerophys(ngrid*nq, zdqssed)
849
850	if(doubleq) then
851	call callsedim2q(ngrid,nlayer, ptimestep,
852	& pplev,zzlev, pt,
853	& pq, pdq, zdqsed, zdqssed,nq)
854	else
855	call callsedim(ngrid,nlayer, ptimestep,
856	& pplev,zzlev, pt,
857	& pq, pdq, zdqsed, zdqssed,nq)
858	end if
859
860
861	DO iq=1, nq
862	DO l=1,nlayer
863	DO ig=1,ngrid
864	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
865	ENDDO
866	ENDDO
867	ENDDO
868	DO iq=1, nq
869	DO ig=1,ngrid
870	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
871	ENDDO
872	ENDDO
873	END IF ! of IF (sedimentation)
874
875	c 7d. Updates
876	c ---------
877
878	DO iq=1, nq
879	DO ig=1,ngrid
880
881	c ---------------------------------
882	c Updating tracer budget on surface
883	c ---------------------------------
884	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
885
886	ENDDO ! (ig)
887	ENDDO ! (iq)
888
889	endif ! of if (tracer)
890
891
892	c-----------------------------------------------------------------------
893	c 8. THERMOSPHERE CALCULATION
894	c-----------------------------------------------------------------------
895
896	if (callthermos) then
897	call thermosphere(pplev,pplay,dist_sol,
898	$ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
899	& pt,pq,pu,pv,pdt,pdq,
900	$ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
901	c do iq=nqchem_min,nq
902	c write(,) 'thermo iq,pq',iq,pq(690,1,iq)
903	c enddo
904
905	DO l=1,nlayer
906	DO ig=1,ngrid
907	dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
908	pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
909	& +zdteuv(ig,l)
910	pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
911	pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
912	DO iq=1, nq
913	pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
914	ENDDO
915	ENDDO
916	ENDDO
917
918
919	endif
920	c-----------------------------------------------------------------------
921	c 9. Surface and sub-surface soil temperature
922	c-----------------------------------------------------------------------
923	c
924	c
925	c 9.1 Increment Surface temperature:
926	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
927
928	DO ig=1,ngrid
929	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
930	ENDDO
931
932	c Prescribe a cold trap at south pole (except at high obliquity !!)
933	c Temperature at the surface is set there to be the temperature
934	c corresponding to equilibrium temperature between phases of CO2
935
936	IF (tracer.AND.water.AND.ngridmx.NE.1) THEN
937	if (caps.and.(obliquit.lt.27.)) then
938	! NB: Updated surface pressure, at grid point 'ngrid', is
939	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
940	tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
941	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
942	endif
943	c -------------------------------------------------------------
944	c Change of surface albedo (set to 0.4) in case of ground frost
945	c everywhere except on the north permanent cap and in regions
946	c covered by dry ice.
947	c ALWAYS PLACE these lines after newcondens !!!
948	c -------------------------------------------------------------
949	do ig=1,ngrid
950
951	c -------------- Version temporaire fit TES 2008 ------------
952	if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then
953	albedo(ig,1)=0.45
954	albedo(ig,2)=0.45
955	endif
956
957	c if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then
958	c if (.not.watercaptag(ig)) then
959	c albedo(ig,1)=0.4
960	c albedo(ig,2)=0.4
961	c endif
962	c endif
963	c -------------- version Francois ---------------
964	c if (co2ice(ig).eq.0.and.
965	c & ((qsurf(ig,nqmx).gt.0.005).or.(watercaptag(ig)))) then
966	c albedo(ig,1)=max(0.4,albedodat(ig))
967	c albedo(ig,2)=albedo(ig,1)
968	c endif
969	c ---------------------------------------------
970	enddo ! of do ig=1,ngrid
971	ENDIF ! of IF (tracer.AND.water.AND.ngridmx.NE.1)
972
973	c
974	c 9.2 Compute soil temperatures and subsurface heat flux:
975	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
976	IF (callsoil) THEN
977	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
978	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
979	ENDIF
980
981	c-----------------------------------------------------------------------
982	c 10. Write output files
983	c ----------------------
984
985	c -------------------------------
986	c Dynamical fields incrementation
987	c -------------------------------
988	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
989	! temperature, zonal and meridional wind
990	DO l=1,nlayer
991	DO ig=1,ngrid
992	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
993	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
994	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
995	ENDDO
996	ENDDO
997
998	! tracers
999	DO iq=1, nq
1000	DO l=1,nlayer
1001	DO ig=1,ngrid
1002	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1003	ENDDO
1004	ENDDO
1005	ENDDO
1006
1007	! surface pressure
1008	DO ig=1,ngrid
1009	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1010	ENDDO
1011
1012	! pressure
1013	DO l=1,nlayer
1014	DO ig=1,ngrid
1015	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1016	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1017	ENDDO
1018	ENDDO
1019
1020	! Density
1021	DO l=1,nlayer
1022	DO ig=1,ngrid
1023	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1024	ENDDO
1025	ENDDO
1026
1027	c Compute surface stress : (NB: z0 is a common in planete.h)
1028	c DO ig=1,ngrid
1029	c cd = (0.4/log(zzlay(ig,1)/z0))**2
1030	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1031	c ENDDO
1032
1033	c Sum of fluxes in solar spectral bands (for output only)
1034	DO ig=1,ngrid
1035	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1036	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1037	ENDDO
1038	c ***** TEST ****************************************************
1039	ztim1 = 999
1040	DO l=1,nlayer
1041	DO ig=1,ngrid
1042	if (pt(ig,l).lt.ztim1) then
1043	ztim1 = pt(ig,l)
1044	igmin = ig
1045	lmin = l
1046	end if
1047	ENDDO
1048	ENDDO
1049	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1050	write(,) 'PHYSIQ: stability WARNING :'
1051	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1052	& 'ig l =', igmin, lmin
1053	end if
1054	c *******************************************************************
1055
1056	c ---------------------
1057	c Outputs to the screen
1058	c ---------------------
1059
1060	IF (lwrite) THEN
1061	PRINT*,'Global diagnostics for the physics'
1062	PRINT,'Variables and their increments x and dx/dt dt'
1063	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1064	WRITE(*,'(2f10.5,2f15.5)')
1065	s tsurf(igout),zdtsurf(igout)*ptimestep,
1066	s pplev(igout,1),pdpsrf(igout)*ptimestep
1067	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1068	WRITE(*,'(i4,6f10.5)') (l,
1069	s pu(igout,l),pdu(igout,l)*ptimestep,
1070	s pv(igout,l),pdv(igout,l)*ptimestep,
1071	s pt(igout,l),pdt(igout,l)*ptimestep,
1072	s l=1,nlayer)
1073	ENDIF ! of IF (lwrite)
1074
1075
1076	IF (ngrid.NE.1) THEN
1077	print,'Ls =',zls180./pi,
1078	& ' tauref(700 Pa,lat=0) =',tauref(ngrid/2),
1079	& ' tau(Viking1) =',tau(ig_vl1,1)
1080
1081	c -------------------------------------------------------------------
1082	c Writing NetCDF file "RESTARTFI" at the end of the run
1083	c -------------------------------------------------------------------
1084	c Note: 'restartfi' is stored just before dynamics are stored
1085	c in 'restart'. Between now and the writting of 'restart',
1086	c there will have been the itau=itau+1 instruction and
1087	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1088	c thus we store for time=time+dtvr
1089
1090	IF(lastcall) THEN
1091	ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1092	write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1093	call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1094	. ptimestep,pday,
1095	. ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
1096	. area,albedodat,inertiedat,zmea,zstd,zsig,
1097	. zgam,zthe)
1098	ENDIF
1099
1100	c -----------------------------------------------------------------
1101	c Saving statistics :
1102	c -----------------------------------------------------------------
1103	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1104	c which can later be used to make the statistic files of the run:
1105	c "stats") only possible in 3D runs !
1106
1107
1108	IF (callstats) THEN
1109
1110
1111	call wstats(ngrid,"ps","Surface pressure","K",2,ps)
1112	call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1113	call wstats(ngrid,"co2ice","CO2 ice cover",
1114	. "kg.m-2",2,co2ice)
1115	c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1116	c . emis)
1117	call wstats(ngrid,"fluxsurf_lw",
1118	. "Thermal IR radiative flux to surface","W.m-2",2,
1119	. fluxsurf_lw)
1120	call wstats(ngrid,"fluxsurf_sw",
1121	. "Solar radiative flux to surface","W.m-2",2,
1122	. fluxsurf_sw_tot)
1123	call wstats(ngrid,"fluxtop_lw",
1124	. "Thermal IR radiative flux to space","W.m-2",2,
1125	. fluxtop_lw)
1126	call wstats(ngrid,"fluxtop_sw",
1127	. "Solar radiative flux to space","W.m-2",2,
1128	. fluxtop_sw_tot)
1129	call wstats(ngrid,"dod","Dust optical depth"," ",2,tau)
1130
1131	call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1132	call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1133	call wstats(ngrid,"v","Meridional (North-South) wind",
1134	. "m.s-1",3,zv)
1135	call wstats(ngrid,"w","Vertical (down-up) wind",
1136	. "m.s-1",3,pw)
1137	call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
1138	call wstats(ngrid,"q2",
1139	. "Boundary layer eddy kinetic energy","m2.s-2",3,q2)
1140
1141	if (tracer) then
1142	if (water) then
1143	vmr=zq(1:ngridmx,1:nlayermx,nqmx)*mugaz/mmol(nqmx)
1144	call wstats(ngrid,"vmr_h2ovapor",
1145	. "H2O vapor volume mixing ratio","mol/mol",
1146	. 3,vmr)
1147	if (iceparty) then
1148	vmr=zq(1:ngridmx,1:nlayermx,nqmx-1)*mugaz/mmol(nqmx-1)
1149	call wstats(ngrid,"vmr_h2oice",
1150	. "H2O ice volume mixing ratio","mol/mol",
1151	. 3,vmr)
1152	endif
1153	endif
1154
1155	if (thermochem.or.photochem) then
1156	do iq=1,nq
1157	if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
1158	. (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
1159	. (noms(iq).eq."h2").or.
1160	. (noms(iq).eq."o3")) then
1161	do l=1,nlayer
1162	do ig=1,ngrid
1163	vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
1164	end do
1165	end do
1166	call wstats(ngrid,"vmr_"//trim(noms(iq)),
1167	. "Volume mixing ratio","mol/mol",3,vmr)
1168	endif
1169	enddo
1170	endif
1171	endif !tracer
1172
1173	IF(lastcall) THEN
1174	write (,) "Writing stats..."
1175	call mkstats(ierr)
1176	ENDIF
1177	ENDIF !if callstats
1178
1179	c (Store EOF for Mars Climate database software)
1180	IF (calleofdump) THEN
1181	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1182	ENDIF
1183
1184
1185	c ----------------------------------------------------------------------
1186	c output in netcdf file "DIAGFI", containing any variable for diagnostic
1187	c (output with period "ecritphy", set in "run.def")
1188	c ----------------------------------------------------------------------
1189	c WRITEDIAGFI can ALSO be called from any other subroutines
1190	c for any variables !!
1191	call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1192	. emis)
1193	call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1194	. tsurf)
1195	call WRITEDIAGFI(ngrid,"ps","surface pressure","K",2,ps)
1196	call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
1197	. co2ice)
1198	c call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1199	c . fluxsurf_lw)
1200	c call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1201	c . fluxsurf_sw_tot)
1202	c call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1203	c . fluxtop_lw)
1204	c call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1205	c . fluxtop_sw_tot)
1206	call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1207	c call WRITEDIAGFI(ngrid,"tau","tau"," ",2,tau)
1208	call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1209	call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1210	call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1211	c call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1212	c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1213	c call WRITEDIAGFI(ngridm,'Teta','T potentielle','K',3,zh)
1214	c call WRITEDIAGFI(ngridm,'Pression','Pression','Pa',3,pplay)
1215	call WRITEDIAGFI(ngrid,"tsoil","soil temperature","K",3,tsoil)
1216
1217
1218	c OUTPUT of tracer mass mixing ratio and surface value :
1219	if (tracer) then
1220	c (for photochemistry, outputs are done in calchim)
1221	do iq=1,nqmx
1222	write(str2(1:2),'(i2.2)') iq
1223	call WRITEDIAGFI(ngridmx,'q'//str2,noms(iq),
1224	& 'kg/kg',3,zq(1,1,iq))
1225	call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq),
1226	& 'kg.m-2',2,qsurf(1,iq))
1227	end do
1228	end if
1229	c ***************************************************
1230
1231	c Outputs for H2O
1232	if (tracer) then
1233	if (activice) then
1234	c call WRITEDIAGFI(ngridmx,'tauice','tau','SI',2,tau(1,2))
1235	c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
1236	c & 'w.m-2',3,zdtsw)
1237	c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
1238	c & 'w.m-2',3,zdtlw)
1239	endif !(activice)
1240
1241	if (water.and..not.photochem) then
1242	iq=nq
1243	c write(str2(1:2),'(i2.2)') iq
1244	c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
1245	c & 'kg.m-2',2,zdqscloud(1,iq))
1246	c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
1247	c & 'kg/kg',3,zdqchim(1,1,iq))
1248	c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
1249	c & 'kg/kg',3,zdqdif(1,1,iq))
1250	c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
1251	c & 'kg/kg',3,zdqadj(1,1,iq))
1252	c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
1253	c & 'kg/kg',3,zdqc(1,1,iq))
1254	endif !(water.and..not.photochem)
1255
1256	c if (iceparty) then
1257	c iq=nq-1
1258	c write(str2(1:2),'(i2.2)') iq
1259	c call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
1260	c & 'kg/kg',3,zq(1,1,iq))
1261	c endif !(iceparty)
1262	endif
1263
1264	c Outputs for dust tracers
1265
1266	if (tracer.and.(dustbin.ne.0)) then
1267	call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
1268	if (doubleq) then
1269	call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
1270	& 'kg.m-2',2,qsurf(1,1))
1271	call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
1272	& 'N.m-2',2,qsurf(1,2))
1273	call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
1274	& 'kg.m-2.s-1',2,zdqsdev(1,1))
1275	call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
1276	& 'kg.m-2.s-1',2,zdqssed(1,1))
1277	do l=1,nlayer
1278	do ig=1, ngrid
1279	reff(ig,l)= ref_r0 *
1280	& (r3n_qpq(ig,l,1)/max(pq(ig,l,2),0.01))*(1./3.)
1281	reff(ig,l)=min(max(reff(ig,l),1.e-10),500.e-6)
1282	end do
1283	end do
1284	call WRITEDIAGFI(ngridmx,'reff','reff','m',3,reff)
1285	else
1286	do iq=1,dustbin
1287	write(str2(1:2),'(i2.2)') iq
1288	call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
1289	& 'kg/kg',3,zq(1,1,iq))
1290	call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
1291	& 'kg.m-2',2,qsurf(1,iq))
1292	end do
1293	endif ! (doubleq)
1294	end if ! (tracer.and.(dustbin.ne.0))
1295
1296
1297	ELSE ! if(ngrid.eq.1)
1298
1299	print,'Ls =',zls180./pi,
1300	& ' tauref(700 Pa) =',tauref,' local tau =',tau(1,1)
1301	c ----------------------------------------------------------------------
1302	c Output in grads file "g1d" (ONLY when using testphys1d)
1303	c (output at every X physical timestep)
1304	c ----------------------------------------------------------------------
1305	c
1306	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
1307	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
1308	CALL writeg1d(ngrid,1,ps,'ps','Pa')
1309
1310	CALL writeg1d(ngrid,nlayer,zt,'T','K')
1311	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
1312	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
1313	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
1314
1315	if(tracer) then
1316	c CALL writeg1d(ngrid,1,tau,'tau','SI')
1317	do iq=1,nq
1318	CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
1319	end do
1320	end if
1321
1322	zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
1323
1324	do l=2,nlayer-1
1325	tmean=zt(1,l)
1326	if(zt(1,l).ne.zt(1,l-1))
1327	& tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
1328	zlocal(l)= zlocal(l-1)
1329	& -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
1330	enddo
1331	zlocal(nlayer)= zlocal(nlayer-1)-
1332	& log(pplay(1,nlayer)/pplay(1,nlayer-1))*
1333	& rnew(1,nlayer)*tmean/g
1334
1335	c if(tracer) then
1336	c do l=2,nlayer
1337	c do iq=1,nq
1338	c hco2(iq)=zq(1,l,iq)/zq(1,l-1,iq)
1339	c hco2(iq)=-(zlocal(l)-zlocal(l-1))/log(hco2(iq))/1000.
1340	c enddo
1341	c write(225,*) l,pt(1,l),(hco2(iq),iq=1,6)
1342	c write(226,*) l,(hco2(iq),iq=7,13)
1343	c enddo
1344	c endif
1345
1346	END IF ! if(ngrid.ne.1)
1347
1348
1349	icount=icount+1
1350	RETURN
1351	END

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

Download in other formats:

Original Format