Context Navigation

physiq.F @ 227

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

MESOSCALE/LMDZ.MARS.new
--> modified to impact last changes

MESOSCALE/LMD_MM_MARS/makemeso
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/call_meso_physiq?.inc
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/call_meso_inifis?.inc
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F
--> modified to get rid of ecri_phys

and make changes related to meso_physiq and meso_inifis

LMDZ.MARS/libf/phymars
--> see LMDZ.MARS/README

15/07/2011 == AS

Modified the mesoscale part so that the previous change by EM does not imply an error in the mesoscale case. More development is needed though to get the "varying z0" capability in the mesoscale model.
Worked on versions of meso_physiq and meso_inifis as close as possible to physiq and inifis for more continuity in the process of impacting changes (and even possibly to reach a common version of physiq and inifis).

The main point is to make the mesoscale significant specific parts

coded into include files in meso_inc so that meso_physiq and meso_inifis looks very close to physiq and inifis.

This is completely transparent for GCM users who does not need the

contents of meso_inc.
Slight cosmetic changes to physiq.f and inifis.F --- some of them e.g. to prepare convergence between meso_physiq and physiq

File size: 57.5 KB

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

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

Download in other formats:

Original Format