Context Navigation

physiq.F @ 173

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

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

Download in other formats:

Original Format