Context Navigation

meso_physiq.F @ 173

Last change on this file since 173 was 173, checked in by acolaitis, 13 years ago
Minor modifications related to thermals
File size: 68.6 KB

Line
1	SUBROUTINE meso_physiq(ngrid,nlayer,nq,
2	$ firstcall,lastcall,
3	$ wday_ini,
4	$ pday,ptime,ptimestep,
5	$ pplev,pplay,pphi,
6	$ pu,pv,pt,pq,pw,
7	$ wtnom,
8	$ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn,
9	$ wtsurf,wtsoil,wemis,wq2,wqsurf,wco2ice,
10	$ wisoil,wdsoil,
11	$ wecritphys,
12	#ifdef MESOSCALE
13	$ output_tab2d, output_tab3d,
14	#endif
15	$ flag_LES)
16
17
18
19	IMPLICIT NONE
20	c=======================================================================
21	c
22	c CAREFUL: THIS IS A VERSION TO BE USED WITH WRF !!!
23	c
24	c ... CHECK THE ****WRF lines
25	c
26	c=======================================================================
27	c
28	c subject:
29	c --------
30	c
31	c Organisation of the physical parametrisations of the LMD
32	c martian atmospheric general circulation model.
33	c
34	c The GCM can be run without or with tracer transport
35	c depending on the value of Logical "tracer" in file "callphys.def"
36	c Tracers may be water vapor, ice OR chemical species OR dust particles
37	c
38	c SEE comments in initracer.F about numbering of tracer species...
39	c
40	c It includes:
41	c
42	c 1. Initialization:
43	c 1.1 First call initializations
44	c 1.2 Initialization for every call to physiq
45	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
46	c 2. Compute radiative transfer tendencies
47	c (longwave and shortwave) for CO2 and aerosols.
48	c 3. Gravity wave and subgrid scale topography drag :
49	c 4. Vertical diffusion (turbulent mixing):
50	c 5. Convective adjustment
51	c 6. Condensation and sublimation of carbon dioxide.
52	c 7. TRACERS :
53	c 7a. water and water ice
54	c 7b. call for photochemistry when tracers are chemical species
55	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
56	c 7d. updates (CO2 pressure variations, surface budget)
57	c 8. Contribution to tendencies due to thermosphere
58	c 9. Surface and sub-surface temperature calculations
59	c 10. Write outputs :
60	c - "startfi", "histfi" (if it's time)
61	c - Saving statistics (if "callstats = .true.")
62	c - Dumping eof (if "calleofdump = .true.")
63	c - Output any needed variables in "diagfi"
64	c 11. Diagnostic: mass conservation of tracers
65	c
66	c author:
67	c -------
68	c Frederic Hourdin 15/10/93
69	c Francois Forget 1994
70	c Christophe Hourdin 02/1997
71	c Subroutine completly rewritten by F.Forget (01/2000)
72	c Introduction of the photochemical module: S. Lebonnois (11/2002)
73	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
74	c Water ice clouds: Franck Montmessin (update 06/2003)
75	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
76	c Nb: See callradite.F for more information.
77	c new WRF version: Aymeric Spiga (01/2009)
78	c
79	c
80	c arguments:
81	c ----------
82	c
83	c input:
84	c ------
85	c ecri period (in dynamical timestep) to write output
86	c ngrid Size of the horizontal grid.
87	c All internal loops are performed on that grid.
88	c nlayer Number of vertical layers.
89	c nq Number of advected fields
90	c firstcall True at the first call
91	c lastcall True at the last call
92	c pday Number of days counted from the North. Spring
93	c equinoxe.
94	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
95	c ptimestep timestep (s)
96	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
97	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
98	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
99	c pu(ngrid,nlayer) u component of the wind (ms-1)
100	c pv(ngrid,nlayer) v component of the wind (ms-1)
101	c pt(ngrid,nlayer) Temperature (K)
102	c pq(ngrid,nlayer,nq) Advected fields
103	c pudyn(ngrid,nlayer) \
104	c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
105	c ptdyn(ngrid,nlayer) / corresponding variables
106	c pqdyn(ngrid,nlayer,nq) /
107	c pw(ngrid,?) vertical velocity
108	c
109	c
110	c ****WRF
111	c day_ini,tsurf,tsoil,emis,q2,qsurf,co2ice are inputs
112	c and locally saved variables
113	c (no need to call phyetat0)
114	c
115	c
116	c output:
117	c -------
118	c
119	c pdu(ngrid,nlayermx) \
120	c pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
121	c pdt(ngrid,nlayermx) / variables due to physical processes.
122	c pdq(ngrid,nlayermx,nqmx) /
123	c pdpsrf(ngrid) /
124	c tracerdyn call tracer in dynamical part of GCM ?
125
126	c
127	c=======================================================================
128	c
129	c 0. Declarations :
130	c ------------------
131
132	#include "dimensions.h"
133	#include "dimphys.h"
134	#include "comgeomfi.h"
135	#include "surfdat.h"
136	#include "comsoil.h" !!! new soil common
137	#include "comdiurn.h"
138	#include "callkeys.h"
139	#include "comcstfi.h"
140	#include "planete.h"
141	#include "comsaison.h"
142	#include "control.h"
143	#include "dimradmars.h"
144	#include "comg1d.h"
145	#include "tracer.h"
146	#include "nlteparams.h"
147
148	#include "chimiedata.h"
149	#include "watercap.h"
150	#include "param.h"
151	#include "param_v3.h"
152	#include "conc.h"
153
154	#include "netcdf.inc"
155
156	!!!!**** SPECIFIC TO MESOSCALE
157	#ifdef MESOSCALE
158	#include "meso_slope.h"
159	#include "wrf_output_2d.h"
160	#include "wrf_output_3d.h"
161	#endif
162
163	#include "advtrac.h" !!! this is necessary for tracers (in dyn3d)
164
165	c Arguments :
166	c -----------
167
168	c inputs:
169	c -------
170	INTEGER ngrid,nlayer,nq
171	REAL ptimestep
172	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
173	REAL pphi(ngridmx,nlayer)
174	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
175	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
176	REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
177	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
178	LOGICAL firstcall,lastcall
179	!!! ****WRF WRF specific to mesoscale
180	INTEGER wday_ini
181	REAL wtsurf(ngridmx) ! input only ay firstcall - output
182	REAL wtsoil(ngridmx,nsoilmx)
183	REAL wisoil(ngridmx,nsoilmx) !! new soil scheme
184	REAL wdsoil(ngridmx,nsoilmx) !! new soil scheme
185	REAL wco2ice(ngridmx)
186	REAL wemis(ngridmx)
187	REAL wqsurf(ngridmx,nqmx)
188	REAL wq2(ngridmx,nlayermx+1)
189	REAL wecritphys
190	#ifdef MESOSCALE
191	REAL output_tab2d(ngridmx,n2d)
192	REAL output_tab3d(ngridmx,nlayer,n3d)
193	#endif
194	REAL sl_ls, sl_lct, sl_lat, sl_tau, sl_alb, sl_the, sl_psi
195	REAL sl_fl0, sl_flu
196	REAL sl_ra, sl_di0
197	REAL sky
198	REAL hfx(ngridmx) !! pour LES avec isfflx!=0
199	REAL ust(ngridmx) !! pour LES avec isfflx!=0
200	LOGICAL flag_LES !! pour LES avec isfflx!=0
201	REAL qsurfice(ngridmx) !! pour diagnostics
202	real alpha,lay1 ! coefficients for building layers
203	integer iloop
204	INTEGER tracerset !!! this corresponds to config%mars
205	!!! ****WRF WRF specific to mesoscale
206	REAL pday
207	REAL ptime
208	logical tracerdyn
209	CHARACTER (len=20) :: wtnom(nqmx) ! tracer name
210
211	c outputs:
212	c --------
213	c physical tendencies
214	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
215	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
216	REAL pdpsrf(ngridmx) ! surface pressure tendency
217
218
219	c Local saved variables:
220	c ----------------------
221	c aerosol (dust or ice) extinction optical depth at reference wavelength
222	c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
223	c aerosol optical properties :
224	REAL aerosol(ngridmx,nlayermx,naerkind)
225
226	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
227	INTEGER icount ! counter of calls to physiq during the run.
228	REAL tsurf(ngridmx) ! Surface temperature (K)
229	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
230	REAL co2ice(ngridmx) ! co2 ice surface layer (kg.m-2)
231	REAL albedo(ngridmx,2) ! Surface albedo in each solar band
232	REAL emis(ngridmx) ! Thermal IR surface emissivity
233	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
234	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
235	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
236	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
237	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
238	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
239	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
240	INTEGER ig_vl1 ! Grid Point near VL1 (for diagnostic)
241
242	c Variables used by the water ice microphysical scheme:
243	REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m)
244	REAL nuice(ngridmx,nlayermx) ! Estimated effective variance
245	! of the size distribution
246	c Albedo of deposited surface ice
247	!!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
248	REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB
249
250	SAVE day_ini, icount
251	SAVE aerosol, tsurf,tsoil
252	SAVE co2ice,albedo,emis, q2
253	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
254	SAVE ig_vl1
255
256	REAL stephan
257	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
258	SAVE stephan
259
260	c Local variables :
261	c -----------------
262
263	REAL CBRT
264	EXTERNAL CBRT
265
266	CHARACTER*80 fichier
267	INTEGER l,ig,ierr,igout,iq,i, tapphys
268
269	REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)
270	REAL fluxsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)
271	REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)
272	REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)
273	REAL tauref(ngridmx) ! Reference column optical depth at 700 Pa
274	! (used if active=F)
275	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
276	REAL zls ! solar longitude (rad)
277	REAL zday ! date (time since Ls=0, in martian days)
278	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
279	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
280	REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
281
282	c Tendancies due to various processes:
283	REAL dqsurf(ngridmx,nqmx)
284	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
285	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
286	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear area
287	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear area
288	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
289	REAL zdtnlte(ngridmx,nlayermx) ! (K/s)
290	REAL zdtsurf(ngridmx) ! (K/s)
291	REAL zdtcloud(ngridmx,nlayermx)
292	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
293	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
294	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
295	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
296	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
297	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
298	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
299	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
300
301	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
302	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
303	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
304	REAL zdqadj(ngridmx,nlayermx,nqmx)
305	REAL zdqc(ngridmx,nlayermx,nqmx)
306	REAL zdqcloud(ngridmx,nlayermx,nqmx)
307	REAL zdqscloud(ngridmx,nqmx)
308	REAL zdqchim(ngridmx,nlayermx,nqmx)
309	REAL zdqschim(ngridmx,nqmx)
310
311	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
312	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
313	REAL zdumolvis(ngridmx,nlayermx)
314	REAL zdvmolvis(ngridmx,nlayermx)
315	real zdqmoldiff(ngridmx,nlayermx,nqmx)
316
317	c Local variable for local intermediate calcul:
318	REAL zflubid(ngridmx)
319	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
320	REAL zdum1(ngridmx,nlayermx)
321	REAL zdum2(ngridmx,nlayermx)
322	REAL ztim1,ztim2,ztim3, z1,z2
323	REAL ztime_fin
324	REAL zdh(ngridmx,nlayermx)
325	INTEGER length
326	PARAMETER (length=100)
327
328	c local variables only used for diagnostic (output in file "diagfi" or "stats")
329	c -----------------------------------------------------------------------------
330	REAL ps(ngridmx), zt(ngridmx,nlayermx)
331	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
332	REAL zq(ngridmx,nlayermx,nqmx)
333	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
334	character*2 str2
335	character*5 str5
336	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
337	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
338	! (particules kg-1)
339	SAVE ccn !! in case iradia != 1
340	real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
341	real qtot1,qtot2 ! total aerosol mass
342	integer igmin, lmin
343	logical tdiag
344
345	real co2col(ngridmx) ! CO2 column
346	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
347	REAL zstress(ngrid), cd
348	real hco2(nqmx),tmean, zlocal(nlayermx)
349	real rho(ngridmx,nlayermx) ! density
350	real vmr(ngridmx,nlayermx) ! volume mixing ratio
351	REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2)
352	REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
353	REAL rave(ngridmx) ! Mean water ice effective radius (m)
354	REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
355	REAL tauTES(ngridmx) ! column optical depth at 825 cm-1
356	REAL Qabsice ! Water ice absorption coefficient
357
358
359	REAL time_phys
360
361	c Variables from thermal
362
363	REAL lmax_th_out(ngrid)
364	REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
365	REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
366	INTEGER lmax_th(ngrid)
367	REAL dtke_th(ngrid,nlayer+1)
368	REAL wmax_th(ngrid),hfmax_th(ngrid)
369	c=======================================================================
370	#ifdef MESOSCALE
371
372	c 1. Initialisation:
373	c -----------------
374
375	c 1.1 Initialisation only at first call
376	c ---------------------------------------
377	IF (firstcall) THEN
378
379	c variables set to 0
380	c ~~~~~~~~~~~~~~~~~~
381	call zerophys(ngridnlayernaerkind,aerosol)
382	call zerophys(ngrid*nlayer,dtrad)
383	call zerophys(ngrid,fluxrad)
384
385	c read startfi
386	c ~~~~~~~~~~~~
387	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
388	c ****WRF
389	c
390	c No need to use startfi.nc
391	c > part of the job of phyetat0 is done in inifis
392	c > remaining initializations are passed here from the WRF variables
393	c > beware, some operations were done by phyetat0 (ex: tracers)
394	c > if any problems, look in phyetat0
395	c
396	tsurf(:)=wtsurf(:)
397	PRINT*,'check: tsurf ',tsurf(1),tsurf(ngridmx)
398	tsoil(:,:)=wtsoil(:,:)
399	PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngridmx,nsoilmx)
400	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
401	!!!new physics
402	inertiedat(:,:)=wisoil(:,:)
403	PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngridmx,nsoilmx)
404	mlayer(0:nsoilmx-1)=wdsoil(1,:)
405	PRINT*,'check: layer ', mlayer
406	!!!!!!!!!!!!!!!!! DONE in soil_setting.F
407	! 1.5 Build layer(); following the same law as mlayer()
408	! Assuming layer distribution follows mid-layer law:
409	! layer(k)=lay1alpha*(k-1)
410	lay1=sqrt(mlayer(0)*mlayer(1))
411	alpha=mlayer(1)/mlayer(0)
412	do iloop=1,nsoilmx
413	layer(iloop)=lay1(alpha*(iloop-1))
414	enddo
415	!!!!!!!!!!!!!!!!! DONE in soil_setting.F
416	tnom(:)=wtnom(:) !! est rempli dans advtrac.h
417	PRINT*,'check: tracernames ', tnom
418	!!!new physics
419	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
420	emis(:)=wemis(:)
421	PRINT*,'check: emis ',emis(1),emis(ngridmx)
422	q2(:,:)=wq2(:,:)
423	PRINT*,'check: q2 ',q2(1,1),q2(ngridmx,nlayermx+1)
424	qsurf(:,:)=wqsurf(:,:)
425	PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngridmx,nqmx)
426	co2ice(:)=wco2ice(:)
427	PRINT*,'check: co2 ',co2ice(1),co2ice(ngridmx)
428	day_ini=wday_ini
429
430	c artificially filling dyn3d/control.h is also required
431	c > iphysiq is put in WRF to be set easily (cf ptimestep)
432	c > day_step is simply deduced:
433	c
434	day_step=daysec/ptimestep
435	PRINT*,'Call to LMD physics:',day_step,' per Martian day'
436	c
437	iphysiq=ptimestep
438	c
439	ecritphy=wecritphys
440	PRINT*,'Write LMD physics each:',ecritphy,' seconds'
441	!!PRINT*,ecri_phys
442	!!PRINT*,float(ecri_phys) ...
443	!!renvoient tous deux des nombres absurdes
444	!!pourtant callkeys.h est inclus ...
445	!!
446	!!donc ecritphys est passe en argument ...
447	PRINT*,'Write LMD physics each:',ecritphy,' seconds'
448	c
449	!DO iq=1, nq
450	! PRINT*, tnom(iq), pq(:,:,iq)
451	!ENDDO
452
453	c
454	c ****WRF
455	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
456
457
458
459	!! Read netcdf initial physical parameters.
460	! CALL phyetat0 ("startfi.nc",0,0,
461	! & nsoilmx,nq,
462	! & day_ini,time_phys,
463	! & tsurf,tsoil,emis,q2,qsurf,co2ice)
464
465	if (pday.ne.day_ini) then
466	write(,) "PHYSIQ: ERROR: bad synchronization between ",
467	& "physics and dynamics"
468	write(,) "dynamics day: ",pday
469	write(,) "physics day: ",day_ini
470	stop
471	endif
472
473	write (,) 'In physiq day_ini =', day_ini
474
475	c Initialize albedo and orbital calculation
476	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
477	CALL surfini(ngrid,co2ice,qsurf,albedo)
478	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
479
480	c initialize soil
481	c ~~~~~~~~~~~~~~~
482	IF (callsoil) THEN
483	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
484	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
485	ELSE
486	PRINT*,
487	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
488	DO ig=1,ngrid
489	capcal(ig)=1.e5
490	fluxgrd(ig)=0.
491	ENDDO
492	ENDIF
493	icount=1
494
495
496	c initialize tracers
497	c ~~~~~~~~~~~~~~~~~~
498	tracerdyn=tracer
499	IF (tracer) THEN
500	CALL initracer(qsurf,co2ice)
501	ENDIF ! end tracer
502
503	!!!!!! WRF WRF WRF MARS MARS
504	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
505	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
506	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
507	!!!!
508	!!!! principe: une option 'caps=T' specifique au mesoscale
509	!!!! ... en vue d'un meso_initracer ????
510	!!!!
511	!!!! depots permanents => albedo TES du PDS
512	!!!! depots saisonniers => alb_surfice (~0.4, cf plus bas)
513	!!!! [!!!! y compris pour les depots saisonniers sur les depots permanents]
514	!!!!
515	!!!! --> todo: il faut garder les depots saisonniers qui viennent
516	!!!! du GCM lorsqu'ils sont consequents
517	!!!!
518	IF ( caps .and. (igcm_h2o_ice .ne. 0) ) THEN
519	PRINT *, 'OVERWRITING watercaptag DEFINITION in INITRACER'
520	PRINT *, 'lat>70 et alb>0.26 => watercaptag=T'
521	!! Perennial H20 north cap defined by watercaptag=true (allows surface to be
522	!! hollowed by sublimation in vdifc).
523	do ig=1,ngridmx
524	qsurf(ig,igcm_h2o_ice)=0. !! on jette les inputs GCM
525	if ( (lati(ig)*180./pi.gt.70.) .and.
526	. (albedodat(ig).ge.0.26) ) then
527	watercaptag(ig)=.true.
528	dryness(ig) = 1.
529	else
530	watercaptag(ig)=.false.
531	dryness(ig) = 1.
532	endif ! (lati, albedodat)
533	end do ! (ngridmx)
534	ELSE ! (caps)
535	print *,'Blork !!!'
536	print *,'caps=T avec water=F ????'
537	ENDIF ! (caps)
538	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
539	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
540	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
541
542
543	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
544	c ****WRF
545	c
546	c nosense in mesoscale modeling
547	c
548	cc Determining gridpoint near Viking Lander 1 (used for diagnostic only)
549	cc ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
550	c
551	c if(ngrid.ne.1) then
552	c latvl1= 22.27
553	c lonvl1= -47.94
554	c ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim +
555	c & int(1.5+(lonvl1+180)*iim/360.)
556	c write(,) 'Viking Lander 1 GCM point: lat,lon',
557	c & lati(ig_vl1)180/pi, long(ig_vl1)180/pi
558	c end if
559	c ****WRF
560	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
561
562	!!!
563	!!! WRF WRF WRF commented for smaller executables
564	!!!
565	!c Initialize thermospheric parameters
566	!c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
567	!
568	! if (callthermos) call param_read
569
570
571	c Initialize R and Cp as constant
572
573	if (.not.callthermos .and. .not.photochem) then
574	do l=1,nlayermx
575	do ig=1,ngridmx
576	rnew(ig,l)=r
577	cpnew(ig,l)=cpp
578	mmean(ig,l)=mugaz
579	enddo
580	enddo
581	endif
582
583	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
584	write(,)"physiq: water_param Surface ice alb:",alb_surfice
585	ENDIF
586
587	ENDIF ! (end of "if firstcall")
588
589	c ---------------------------------------------------
590	c 1.2 Initializations done at every physical timestep:
591	c ---------------------------------------------------
592	c
593	IF (ngrid.NE.ngridmx) THEN
594	PRINT*,'STOP in PHYSIQ'
595	PRINT*,'Probleme de dimensions :'
596	PRINT*,'ngrid = ',ngrid
597	PRINT*,'ngridmx = ',ngridmx
598	STOP
599	ENDIF
600
601	c Initialize various variables
602	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
603	call zerophys(ngrid*nlayer, pdv)
604	call zerophys(ngrid*nlayer, pdu)
605	call zerophys(ngrid*nlayer, pdt)
606	call zerophys(ngridnlayernq, pdq)
607	call zerophys(ngrid, pdpsrf)
608	call zerophys(ngrid, zflubid)
609	call zerophys(ngrid, zdtsurf)
610	call zerophys(ngrid*nq, dqsurf)
611	igout=ngrid/2+1
612
613
614	zday=pday+ptime ! compute time, in sols (and fraction thereof)
615
616	c Compute Solar Longitude (Ls) :
617	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
618	if (season) then
619	call solarlong(zday,zls)
620	else
621	call solarlong(float(day_ini),zls)
622	end if
623
624	c Compute geopotential at interlayers
625	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
626	c ponderation des altitudes au niveau des couches en dp/p
627
628	DO l=1,nlayer
629	DO ig=1,ngrid
630	zzlay(ig,l)=pphi(ig,l)/g
631	ENDDO
632	ENDDO
633	DO ig=1,ngrid
634	zzlev(ig,1)=0.
635	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
636	ENDDO
637	DO l=2,nlayer
638	DO ig=1,ngrid
639	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
640	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
641	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
642	ENDDO
643	ENDDO
644
645
646	! Potential temperature calculation not the same in physiq and dynamic
647
648	c Compute potential temperature
649	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
650	DO l=1,nlayer
651	DO ig=1,ngrid
652	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
653	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
654	ENDDO
655	ENDDO
656
657	!!!
658	!!! WRF WRF WRF commented for smaller executables
659	!!!
660	!c-----------------------------------------------------------------------
661	!c 1.2.5 Compute mean mass, cp, and R
662	!c --------------------------------
663	!
664	! if(photochem.or.callthermos) then
665	! call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
666	! endif
667
668	c-----------------------------------------------------------------------
669	c 2. Compute radiative tendencies :
670	c------------------------------------
671
672
673	IF (callrad) THEN
674	IF( MOD(icount-1,iradia).EQ.0) THEN
675
676	write (,) 'call radiative transfer'
677
678	c Local Solar zenith angle
679	c ~~~~~~~~~~~~~~~~~~~~~~~~
680	CALL orbite(zls,dist_sol,declin)
681
682	IF(diurnal) THEN
683	ztim1=SIN(declin)
684	ztim2=COS(declin)COS(2.pi*(zday-.5))
685	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
686
687	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
688	s ztim1,ztim2,ztim3, mu0,fract)
689
690	ELSE
691	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
692	ENDIF
693
694	c NLTE cooling from CO2 emission
695	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
696
697	IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
698
699	c Find number of layers for LTE radiation calculations
700	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
701	& CALL nlthermeq(ngrid,nlayer,pplev,pplay)
702
703	c Note: Dustopacity.F has been transferred to callradite.F
704
705	c Call main radiative transfer scheme
706	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
707	c Transfer through CO2 (except NIR CO2 absorption)
708	c and aerosols (dust and water ice)
709
710	c Radiative transfer
711	c ------------------
712	cc
713	cc **WRF: desormais dust_opacity est dans callradite -- modifications
714	cc nveaux arguments: tauref,tau,aerosol,rice,nuice
715	cc
716	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
717	$ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
718	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
719	& tauref,tau,aerosol,ccn,rdust,rice,nuice)
720
721	c write(,) icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
722	c $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
723	c $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
724	c & tauref,tau,aerosol,rice,nuice
725	c write(,) fluxsurf_lw
726
727	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
728	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
729	ccccc
730	ccccc PARAM SLOPE : Insolation (direct + scattered)
731	ccccc
732	DO ig=1,ngrid
733	sl_the = theta_sl(ig)
734	IF (sl_the .ne. 0.) THEN
735	ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
736	DO l=1,2
737	sl_lct = ptime24. + 180.long(ig)/pi/15.
738	sl_ra = pi*(1.0-sl_lct/12.)
739	sl_lat = 180.*lati(ig)/pi
740	sl_tau = tau(ig,1)
741	sl_alb = albedo(ig,l)
742	sl_psi = psi_sl(ig)
743	sl_fl0 = fluxsurf_sw(ig,l)
744	sl_di0 = 0.
745	if (mu0(ig) .gt. 0.) then
746	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
747	sl_di0 = sl_di0*1370./dist_sol/dist_sol
748	sl_di0 = sl_di0/ztim1
749	sl_di0 = fluxsurf_sw(ig,l)*sl_di0
750	endif
751	! sait-on jamais (a cause des arrondis)
752	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
753	!!!!!!!!!!!!!!!!!!!!!!!!!!
754	CALL meso_param_slope( mu0(ig), declin, sl_ra, sl_lat,
755	& sl_tau, sl_alb,
756	& sl_the, sl_psi, sl_di0, sl_fl0, sl_flu)
757	!!!!!!!!!!!!!!!!!!!!!!!!!!
758	fluxsurf_sw(ig,l) = sl_flu
759	!! sl_ls = 180.*zls/pi
760	!! sl_lct = ptime24. + 180.long(ig)/pi/15.
761	!! sl_lat = 180.*lati(ig)/pi
762	!! sl_tau = tau(ig,1)
763	!! sl_alb = albedo(ig,l)
764	!! sl_the = theta_sl(ig)
765	!! sl_psi = psi_sl(ig)
766	!! sl_fl0 = fluxsurf_sw(ig,l)
767	!! CALL param_slope_full(sl_ls, sl_lct, sl_lat,
768	!! & sl_tau, sl_alb,
769	!! & sl_the, sl_psi, sl_fl0, sl_flu)
770	ENDDO
771	!!! compute correction on IR flux as well
772	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
773	fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
774	ENDIF
775	ENDDO
776	ccccc
777	ccccc PARAM SLOPE
778	ccccc
779	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
780	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
781
782
783	c CO2 near infrared absorption
784	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
785	call zerophys(ngrid*nlayer,zdtnirco2)
786	if (callnirco2) then
787	call nirco2abs (ngrid,nlayer,pplay,dist_sol,
788	. mu0,fract,declin, zdtnirco2)
789	endif
790
791	c Radiative flux from the sky absorbed by the surface (W.m-2)
792	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
793	DO ig=1,ngrid
794	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
795	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
796	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
797
798	!print*,'RAD ', fluxrad_sky(ig)
799	!print,'LW ', emis(ig)fluxsurf_lw(ig)
800	!print,'SW ', fluxsurf_sw(ig,2)(1.-albedo(ig,2))
801
802	ENDDO
803
804
805	c Net atmospheric radiative heating rate (K.s-1)
806	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
807	IF(callnlte) THEN
808	CALL blendrad(ngrid, nlayer, pplay,
809	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
810	ELSE
811	DO l=1,nlayer
812	DO ig=1,ngrid
813	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
814	& +zdtnirco2(ig,l)
815	ENDDO
816	ENDDO
817	ENDIF
818
819
820
821	ENDIF ! of if(mod(icount-1,iradia).eq.0)
822
823	c Transformation of the radiative tendencies:
824	c -------------------------------------------
825
826	c Net radiative surface flux (W.m-2)
827	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
828	c
829	DO ig=1,ngrid
830	zplanck(ig)=tsurf(ig)*tsurf(ig)
831	zplanck(ig)=emis(ig)*
832	$ stephanzplanck(ig)zplanck(ig)
833	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
834	cccc
835	cccc param slope
836	cccc
837	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
838	fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
839	cccc
840	cccc
841	cccc
842	ENDDO
843
844
845	DO l=1,nlayer
846	DO ig=1,ngrid
847	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
848	ENDDO
849	ENDDO
850
851	ENDIF ! of IF (callrad)
852
853	!!!
854	!!! WRF WRF WRF commented for smaller executables
855	!!!
856	!c-----------------------------------------------------------------------
857	!c 3. Gravity wave and subgrid scale topography drag :
858	!c -------------------------------------------------
859	!
860	!
861	! IF(calllott)THEN
862	!
863	! CALL calldrag_noro(ngrid,nlayer,ptimestep,
864	! & pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
865	!
866	! DO l=1,nlayer
867	! DO ig=1,ngrid
868	! pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
869	! pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
870	! pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
871	! ENDDO
872	! ENDDO
873	! ENDIF
874
875	c-----------------------------------------------------------------------
876	c 4. Vertical diffusion (turbulent mixing):
877	c -----------------------------------------
878	c
879	IF (calldifv) THEN
880
881
882	DO ig=1,ngrid
883	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
884	!write (,), fluxrad(ig), fluxgrd(ig), zflubid(ig)
885	ENDDO
886
887	CALL zerophys(ngrid*nlayer,zdum1)
888	CALL zerophys(ngrid*nlayer,zdum2)
889	do l=1,nlayer
890	do ig=1,ngrid
891	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
892	enddo
893	enddo
894
895	c Calling vdif (Martian version WITH CO2 condensation)
896	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
897	$ ptimestep,capcal,lwrite,
898	$ pplay,pplev,zzlay,zzlev,z0,
899	$ pu,pv,zh,pq,tsurf,emis,qsurf,
900	$ zdum1,zdum2,zdh,pdq,zflubid,
901	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
902	& zdqdif,zdqsdif)
903
904	DO ig=1,ngrid
905	!! sensible heat flux in W/m2
906	hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
907	!! u star in similarity theory in m/s
908	ust(ig) = 0.4
909	. * sqrt( pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1) )
910	. / log( 1.E+0 + zzlay(ig,1)/z0 )
911	ENDDO
912
913	! write (,) 'PHYS HFX cp zdts', hfx(100), zflubid(100),
914	! . capcal(100),
915	! . zdtsdif(100)
916	! write (,) 'PHYS UST', ust(100)
917
918	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
919	!!! LES LES
920	IF (flag_LES) THEN
921
922	write (,) '************************************************'
923	write (,) '** LES mode: the difv part is only used to'
924	write (,) '** provide HFX and UST to the dynamics'
925	write (,) '** NB: - dudif, dvdif, dhdif, dqdif are set to 0'
926	write (,) '** - tsurf is updated'
927	write (,) '************************************************'
928
929	!!!
930	!!! WRF WRF LES LES : en fait le subgrid scale n'etait pas mis a zero !!
931	!!!
932	DO ig=1,ngrid
933	! !! sensible heat flux in W/m2
934	! hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
935	! !! u star in similarity theory in m/s
936	! ust(ig) = 0.4
937	! . * sqrt( pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1) )
938	! . / log( 1.E+0 + zzlay(ig,1)/z0 )
939	!
940	DO l=1,nlayer
941	zdvdif(ig,l) = 0.
942	zdudif(ig,l) = 0.
943	zdhdif(ig,l) = 0.
944	DO iq=1, nq
945	zdqdif(ig,l,iq) = 0.
946	zdqsdif(ig,iq) = 0. !! sortir de la boucle
947	ENDDO
948	ENDDO
949	!
950	ENDDO
951	!write (,) 'RAD ',fluxrad(igout)
952	!write (,) 'GRD ',fluxgrd(igout)
953	!write (,) 'dTs/dt ',capcal(igout)*zdtsurf(igout)
954	!write (,) 'HFX ', hfx(igout)
955	!write (,) 'UST ', ust(igout)
956	ENDIF
957	!!! LES LES
958	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
959
960	DO l=1,nlayer
961	DO ig=1,ngrid
962	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
963	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
964	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
965
966	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
967
968	ENDDO
969	ENDDO
970
971	DO ig=1,ngrid
972	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
973	ENDDO
974
975	if (tracer) then
976	DO iq=1, nq
977	DO l=1,nlayer
978	DO ig=1,ngrid
979	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
980	ENDDO
981	ENDDO
982	ENDDO
983	DO iq=1, nq
984	DO ig=1,ngrid
985	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
986	ENDDO
987	ENDDO
988	end if ! of if (tracer)
989
990	ELSE
991	DO ig=1,ngrid
992	zdtsurf(ig)=zdtsurf(ig)+
993	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
994	ENDDO
995	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
996	IF (flag_LES) THEN
997	write(,) 'LES mode !'
998	write(,) 'Please set calldifv to T in callphys.def'
999	STOP
1000	ENDIF
1001	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1002	ENDIF ! of IF (calldifv)
1003
1004	c-----------------------------------------------------------------------
1005	c TEST. Thermals :
1006	c HIGHLY EXPERIMENTAL, BEWARE !!
1007	c -----------------------------
1008
1009	if(calltherm) then
1010
1011	call calltherm_interface(ngrid,nlayer,firstcall,
1012	$ long,lati,zzlev,zzlay,
1013	$ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
1014	$ pplay,pplev,pphi,nq,zpopsk,
1015	$ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,dtke_th,hfmax_th,wmax_th)
1016
1017	DO l=1,nlayer
1018	DO ig=1,ngrid
1019	pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
1020	pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
1021	pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
1022	q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
1023	ENDDO
1024	ENDDO
1025
1026	DO ig=1,ngrid
1027	q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
1028	ENDDO
1029
1030	if (tracer) then
1031	DO iq=1,nq
1032	DO l=1,nlayer
1033	DO ig=1,ngrid
1034	pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
1035	ENDDO
1036	ENDDO
1037	ENDDO
1038	endif
1039
1040	else !of if calltherm
1041	lmax_th(:)=0
1042	end if
1043
1044	c-----------------------------------------------------------------------
1045	c 5. Dry convective adjustment:
1046	c -----------------------------
1047
1048	IF(calladj) THEN
1049
1050	DO l=1,nlayer
1051	DO ig=1,ngrid
1052	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
1053	ENDDO
1054	ENDDO
1055	CALL zerophys(ngrid*nlayer,zduadj)
1056	CALL zerophys(ngrid*nlayer,zdvadj)
1057	CALL zerophys(ngrid*nlayer,zdhadj)
1058
1059	CALL convadj(ngrid,nlayer,nq,ptimestep,
1060	$ pplay,pplev,zpopsk,lmax_th,
1061	$ pu,pv,zh,pq,
1062	$ pdu,pdv,zdh,pdq,
1063	$ zduadj,zdvadj,zdhadj,
1064	$ zdqadj)
1065
1066
1067	DO l=1,nlayer
1068	DO ig=1,ngrid
1069	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
1070	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
1071	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
1072
1073	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
1074	ENDDO
1075	ENDDO
1076
1077	if(tracer) then
1078	DO iq=1, nq
1079	DO l=1,nlayer
1080	DO ig=1,ngrid
1081	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
1082	ENDDO
1083	ENDDO
1084	ENDDO
1085	end if
1086	ENDIF ! of IF(calladj)
1087
1088	c-----------------------------------------------------------------------
1089	c 6. Carbon dioxide condensation-sublimation:
1090	c -------------------------------------------
1091
1092	IF (callcond) THEN
1093	CALL newcondens(ngrid,nlayer,nq,ptimestep,
1094	$ capcal,pplay,pplev,tsurf,pt,
1095	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
1096	$ co2ice,albedo,emis,
1097	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
1098	$ fluxsurf_sw,zls)
1099
1100	DO l=1,nlayer
1101	DO ig=1,ngrid
1102	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
1103	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
1104	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
1105	ENDDO
1106	ENDDO
1107	DO ig=1,ngrid
1108	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
1109	ENDDO
1110
1111	IF (tracer) THEN
1112	DO iq=1, nq
1113	DO l=1,nlayer
1114	DO ig=1,ngrid
1115	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1116	ENDDO
1117	ENDDO
1118	ENDDO
1119	ENDIF ! of IF (tracer)
1120
1121	ENDIF ! of IF (callcond)
1122
1123	c-----------------------------------------------------------------------
1124	c 7. Specific parameterizations for tracers
1125	c: -----------------------------------------
1126
1127	if (tracer) then
1128
1129	c 7a. Water and ice
1130	c ---------------
1131
1132	c ---------------------------------------
1133	c Water ice condensation in the atmosphere
1134	c ----------------------------------------
1135	IF (water) THEN
1136
1137	call watercloud(ngrid,nlayer,ptimestep,
1138	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
1139	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
1140	& nq,naerkind,tau,
1141	& ccn,rdust,rice,nuice)
1142	if (activice) then
1143	c Temperature variation due to latent heat release
1144	DO l=1,nlayer
1145	DO ig=1,ngrid
1146	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
1147	ENDDO
1148	ENDDO
1149	endif
1150
1151	! increment water vapour and ice atmospheric tracers tendencies
1152	IF (water) THEN
1153	DO l=1,nlayer
1154	DO ig=1,ngrid
1155	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
1156	& zdqcloud(ig,l,igcm_h2o_vap)
1157	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
1158	& zdqcloud(ig,l,igcm_h2o_ice)
1159	ENDDO
1160	ENDDO
1161	ENDIF ! of IF (water) THEN
1162	! Increment water ice surface tracer tendency
1163	DO ig=1,ngrid
1164	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
1165	& zdqscloud(ig,igcm_h2o_ice)
1166	ENDDO
1167
1168	END IF ! of IF (water)
1169
1170	c 7b. Chemical species
1171	c ------------------
1172
1173	!!!
1174	!!! WRF WRF WRF commented for smaller executables
1175	!!!
1176	!c --------------
1177	!c photochemistry :
1178	!c --------------
1179	! IF (photochem .or. thermochem) then
1180	! call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
1181	! & zzlay,zday,pq,pdq,rice,
1182	! & zdqchim,zdqschim,zdqcloud,zdqscloud)
1183	!!NB: Photochemistry includes condensation of H2O2
1184	!
1185	! ! increment values of tracers:
1186	! DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1187	! ! tracers is zero anyways
1188	! DO l=1,nlayer
1189	! DO ig=1,ngrid
1190	! pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
1191	! ENDDO
1192	! ENDDO
1193	! ENDDO ! of DO iq=1,nq
1194	! ! add condensation tendency for H2O2
1195	! if (igcm_h2o2.ne.0) then
1196	! DO l=1,nlayer
1197	! DO ig=1,ngrid
1198	! pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
1199	! & +zdqcloud(ig,l,igcm_h2o2)
1200	! ENDDO
1201	! ENDDO
1202	! endif
1203	!
1204	! ! increment surface values of tracers:
1205	! DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1206	! ! tracers is zero anyways
1207	! DO ig=1,ngrid
1208	! dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
1209	! ENDDO
1210	! ENDDO ! of DO iq=1,nq
1211	! ! add condensation tendency for H2O2
1212	! if (igcm_h2o2.ne.0) then
1213	! DO ig=1,ngrid
1214	! dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
1215	! & +zdqscloud(ig,igcm_h2o2)
1216	! ENDDO
1217	! endif
1218	!
1219	! END IF ! of IF (photochem.or.thermochem)
1220
1221	c 7c. Aerosol particles
1222	c -------------------
1223
1224	c ----------
1225	c Dust devil :
1226	c ----------
1227	IF(callddevil) then
1228	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
1229	& zdqdev,zdqsdev)
1230
1231	if (dustbin.ge.1) then
1232	do iq=1,nq
1233	DO l=1,nlayer
1234	DO ig=1,ngrid
1235	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
1236	ENDDO
1237	ENDDO
1238	enddo
1239	do iq=1,nq
1240	DO ig=1,ngrid
1241	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
1242	ENDDO
1243	enddo
1244	endif ! of if (dustbin.ge.1)
1245
1246	END IF ! of IF (callddevil)
1247
1248	c -------------
1249	c Sedimentation : acts also on water ice
1250	c -------------
1251	IF (sedimentation) THEN
1252	!call zerophys(ngridnlayernq, zdqsed)
1253	zdqsed(1:ngrid,1:nlayer,1:nq)=0
1254	!call zerophys(ngrid*nq, zdqssed)
1255	zdqssed(1:ngrid,1:nq)=0
1256
1257	call callsedim(ngrid,nlayer, ptimestep,
1258	& pplev,zzlev, pt, rdust, rice,
1259	& pq, pdq, zdqsed, zdqssed,nq)
1260
1261	DO iq=1, nq
1262	DO l=1,nlayer
1263	DO ig=1,ngrid
1264	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1265	ENDDO
1266	ENDDO
1267	ENDDO
1268	DO iq=1, nq
1269	DO ig=1,ngrid
1270	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1271	ENDDO
1272	ENDDO
1273	END IF ! of IF (sedimentation)
1274
1275	c 7d. Updates
1276	c ---------
1277
1278	DO iq=1, nq
1279	DO ig=1,ngrid
1280
1281	c ---------------------------------
1282	c Updating tracer budget on surface
1283	c ---------------------------------
1284	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1285
1286	ENDDO ! (ig)
1287	ENDDO ! (iq)
1288
1289	endif ! of if (tracer)
1290
1291	!!!
1292	!!! WRF WRF WRF commented for smaller executables
1293	!!!
1294	!c-----------------------------------------------------------------------
1295	!c 8. THERMOSPHERE CALCULATION
1296	!c-----------------------------------------------------------------------
1297	!
1298	! if (callthermos) then
1299	! call thermosphere(pplev,pplay,dist_sol,
1300	! $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1301	! & pt,pq,pu,pv,pdt,pdq,
1302	! $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1303	!
1304	! DO l=1,nlayer
1305	! DO ig=1,ngrid
1306	! dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1307	! pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1308	! & +zdteuv(ig,l)
1309	! pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1310	! pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1311	! DO iq=1, nq
1312	! pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1313	! ENDDO
1314	! ENDDO
1315	! ENDDO
1316	!
1317	! endif ! of if (callthermos)
1318
1319	c-----------------------------------------------------------------------
1320	c 9. Surface and sub-surface soil temperature
1321	c-----------------------------------------------------------------------
1322	c
1323	c
1324	c 9.1 Increment Surface temperature:
1325	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1326
1327	DO ig=1,ngrid
1328	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1329	ENDDO
1330
1331	c Prescribe a cold trap at south pole (except at high obliquity !!)
1332	c Temperature at the surface is set there to be the temperature
1333	c corresponding to equilibrium temperature between phases of CO2
1334
1335	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
1336	! if (caps.and.(obliquit.lt.27.)) then
1337	! ! NB: Updated surface pressure, at grid point 'ngrid', is
1338	! ! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1339	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1340	! & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1341	! endif
1342	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1343	!!!!! note WRF MESOSCALE AYMERIC -- mot cle "caps"
1344	!!!!! watercaptag n'est plus utilise que dans vdifc
1345	!!!!! ... pour que la sublimation ne soit pas stoppee
1346	!!!!! ... dans la calotte permanente nord si qsurf=0
1347	!!!!! on desire garder cet effet regle par caps=T
1348	!!!!! on a donc commente "if (caps.and.(obliquit.lt.27.))" ci-dessus
1349	!!!!! --- remplacer ces lignes par qqch de plus approprie
1350	!!!!! si on s attaque a la calotte polaire sud
1351	!!!!! pas d'autre occurrence majeure du mot-cle "caps"
1352	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1353
1354	c -------------------------------------------------------------
1355	c Change of surface albedo (set to 0.4) in case of ground frost
1356	c everywhere except on the north permanent cap and in regions
1357	c covered by dry ice.
1358	c ALWAYS PLACE these lines after newcondens !!!
1359	c -------------------------------------------------------------
1360	c **WRF : OK avec le mesoscale, pas d'indices bizarres au pole
1361	do ig=1,ngrid
1362	if ((co2ice(ig).eq.0).and.
1363	& (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
1364	albedo(ig,1) = alb_surfice
1365	albedo(ig,2) = alb_surfice
1366	endif
1367	enddo ! of do ig=1,ngrid
1368	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
1369
1370	c
1371	c 9.2 Compute soil temperatures and subsurface heat flux:
1372	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1373	IF (callsoil) THEN
1374	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1375	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
1376	ENDIF
1377
1378	c-----------------------------------------------------------------------
1379	c 10. Write output files
1380	c ----------------------
1381
1382	c -------------------------------
1383	c Dynamical fields incrementation
1384	c -------------------------------
1385	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1386	! temperature, zonal and meridional wind
1387	DO l=1,nlayer
1388	DO ig=1,ngrid
1389	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1390	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1391	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1392	ENDDO
1393	ENDDO
1394
1395	! tracers
1396	DO iq=1, nq
1397	DO l=1,nlayer
1398	DO ig=1,ngrid
1399	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1400	ENDDO
1401	ENDDO
1402	ENDDO
1403
1404	! surface pressure
1405	DO ig=1,ngrid
1406	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1407	ENDDO
1408
1409	! pressure
1410	DO l=1,nlayer
1411	DO ig=1,ngrid
1412	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1413	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1414	ENDDO
1415	ENDDO
1416
1417	! Density
1418	DO l=1,nlayer
1419	DO ig=1,ngrid
1420	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1421	ENDDO
1422	ENDDO
1423
1424	c Compute surface stress : (NB: z0 is a common in planete.h)
1425	c DO ig=1,ngrid
1426	c cd = (0.4/log(zzlay(ig,1)/z0))**2
1427	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1428	c ENDDO
1429
1430	c Sum of fluxes in solar spectral bands (for output only)
1431	DO ig=1,ngrid
1432	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1433	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1434	ENDDO
1435	c ***** TEST ****************************************************
1436	ztim1 = 999
1437	DO l=1,nlayer
1438	DO ig=1,ngrid
1439	if (pt(ig,l).lt.ztim1) then
1440	ztim1 = pt(ig,l)
1441	igmin = ig
1442	lmin = l
1443	end if
1444	ENDDO
1445	ENDDO
1446	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1447	write(,) 'PHYSIQ: stability WARNING :'
1448	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1449	& 'ig l =', igmin, lmin
1450	end if
1451	c *******************************************************************
1452
1453	c ---------------------
1454	c Outputs to the screen
1455	c ---------------------
1456
1457	IF (lwrite) THEN
1458	PRINT*,'Global diagnostics for the physics'
1459	PRINT,'Variables and their increments x and dx/dt dt'
1460	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1461	WRITE(*,'(2f10.5,2f15.5)')
1462	s tsurf(igout),zdtsurf(igout)*ptimestep,
1463	s pplev(igout,1),pdpsrf(igout)*ptimestep
1464	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1465	WRITE(*,'(i4,6f10.5)') (l,
1466	s pu(igout,l),pdu(igout,l)*ptimestep,
1467	s pv(igout,l),pdv(igout,l)*ptimestep,
1468	s pt(igout,l),pdt(igout,l)*ptimestep,
1469	s l=1,nlayer)
1470	ENDIF ! of IF (lwrite)
1471
1472	IF (ngrid.NE.1) THEN
1473	print,'Ls =',zls180./pi,
1474	& ' tauref(700 Pa,lat=0) =',tauref(ngrid/2)!,
1475	! & ' tau(Viking1) =',tau(ig_vl1,1)
1476
1477
1478	c -------------------------------------------------------------------
1479	c Writing NetCDF file "RESTARTFI" at the end of the run
1480	c -------------------------------------------------------------------
1481	c Note: 'restartfi' is stored just before dynamics are stored
1482	c in 'restart'. Between now and the writting of 'restart',
1483	c there will have been the itau=itau+1 instruction and
1484	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1485	c thus we store for time=time+dtvr
1486
1487
1488	!!!
1489	!!! WRF WRF WRF WRF
1490	!!!
1491	! IF(lastcall) THEN
1492	! ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1493	! write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1494	! call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1495	! . ptimestep,pday,
1496	! . ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
1497	! . area,albedodat,inertiedat,zmea,zstd,zsig,
1498	! . zgam,zthe)
1499	! ENDIF
1500
1501
1502
1503	c -------------------------------------------------------------------
1504	c Calculation of diagnostic variables written in both stats and
1505	c diagfi files
1506	c -------------------------------------------------------------------
1507
1508	if (tracer) then
1509	if (water) then
1510
1511	call zerophys(ngrid,mtot)
1512	call zerophys(ngrid,icetot)
1513	call zerophys(ngrid,rave)
1514	call zerophys(ngrid,tauTES)
1515	do ig=1,ngrid
1516	do l=1,nlayermx
1517	mtot(ig) = mtot(ig) +
1518	& zq(ig,l,igcm_h2o_vap) *
1519	& (pplev(ig,l) - pplev(ig,l+1)) / g
1520	icetot(ig) = icetot(ig) +
1521	& zq(ig,l,igcm_h2o_ice) *
1522	& (pplev(ig,l) - pplev(ig,l+1)) / g
1523	rave(ig) = rave(ig) +
1524	& zq(ig,l,igcm_h2o_ice) *
1525	& (pplev(ig,l) - pplev(ig,l+1)) / g *
1526	& rice(ig,l) * (1.+nuice_ref)
1527	c Computing abs optical depth at 825 cm-1 in each
1528	c layer to simulate NEW TES retrieval
1529	Qabsice = min(
1530	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
1531	& )
1532	opTES(ig,l)= 0.75 * Qabsice *
1533	& zq(ig,l,igcm_h2o_ice) *
1534	& (pplev(ig,l) - pplev(ig,l+1)) / g
1535	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1536	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1537	enddo
1538	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
1539	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1540	enddo
1541
1542	endif ! of if (water)
1543	endif ! of if (tracer)
1544
1545	c -----------------------------------------------------------------
1546	c WSTATS: Saving statistics
1547	c -----------------------------------------------------------------
1548	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1549	c which can later be used to make the statistic files of the run:
1550	c "stats") only possible in 3D runs !
1551
1552	IF (callstats) THEN
1553
1554	write(,) 'callstats'
1555
1556	! call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1557	! call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1558	! call wstats(ngrid,"co2ice","CO2 ice cover",
1559	! & "kg.m-2",2,co2ice)
1560	! call wstats(ngrid,"fluxsurf_lw",
1561	! & "Thermal IR radiative flux to surface","W.m-2",2,
1562	! & fluxsurf_lw)
1563	! call wstats(ngrid,"fluxsurf_sw",
1564	! & "Solar radiative flux to surface","W.m-2",2,
1565	! & fluxsurf_sw_tot)
1566	! call wstats(ngrid,"fluxtop_lw",
1567	! & "Thermal IR radiative flux to space","W.m-2",2,
1568	! & fluxtop_lw)
1569	! call wstats(ngrid,"fluxtop_sw",
1570	! & "Solar radiative flux to space","W.m-2",2,
1571	! & fluxtop_sw_tot)
1572	! call wstats(ngrid,"taudustvis",
1573	! & "Dust optical depth"," ",2,tau(1,1))
1574	! call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1575	! call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1576	! call wstats(ngrid,"v","Meridional (North-South) wind",
1577	! & "m.s-1",3,zv)
1578	!c call wstats(ngrid,"w","Vertical (down-up) wind",
1579	!c & "m.s-1",3,pw)
1580	! call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
1581	!c call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
1582	!c call wstats(ngrid,"q2",
1583	!c & "Boundary layer eddy kinetic energy",
1584	!c & "m2.s-2",3,q2)
1585	!c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1586	!c & emis)
1587	!c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
1588	!c & 2,zstress)
1589	!
1590	! if (tracer) then
1591	! if (water) then
1592	! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1593	! & *mugaz/mmol(igcm_h2o_vap)
1594	! call wstats(ngrid,"vmr_h2ovapor",
1595	! & "H2O vapor volume mixing ratio","mol/mol",
1596	! & 3,vmr)
1597	! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1598	! & *mugaz/mmol(igcm_h2o_ice)
1599	! call wstats(ngrid,"vmr_h2oice",
1600	! & "H2O ice volume mixing ratio","mol/mol",
1601	! & 3,vmr)
1602	!
1603	! call wstats(ngrid,"mtot",
1604	! & "total mass of water vapor","kg/m2",
1605	! & 2,mtot)
1606	! call wstats(ngrid,"icetot",
1607	! & "total mass of water ice","kg/m2",
1608	! & 2,icetot)
1609	!c If activice is true, tauTES is computed in aeropacity.F;
1610	! if (.not.activice) then
1611	! call wstats(ngrid,"tauTES",
1612	! & "tau abs 825 cm-1","",
1613	! & 2,tauTES)
1614	! endif ! of if (activice)
1615	!
1616	! endif ! of if (water)
1617	!
1618	! if (thermochem.or.photochem) then
1619	! do iq=1,nq
1620	! if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
1621	! . (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
1622	! . (noms(iq).eq."h2").or.
1623	! . (noms(iq).eq."o3")) then
1624	! do l=1,nlayer
1625	! do ig=1,ngrid
1626	! vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
1627	! end do
1628	! end do
1629	! call wstats(ngrid,"vmr_"//trim(noms(iq)),
1630	! . "Volume mixing ratio","mol/mol",3,vmr)
1631	! endif
1632	! enddo
1633	! endif ! of if (thermochem.or.photochem)
1634	!
1635	! endif ! of if (tracer)
1636	!
1637	! IF(lastcall) THEN
1638	! write (,) "Writing stats..."
1639	! call mkstats(ierr)
1640	! ENDIF
1641
1642	ENDIF !if callstats
1643
1644	c (Store EOF for Mars Climate database software)
1645	IF (calleofdump) THEN
1646	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1647	ENDIF
1648
1649	ccc************** WRF OUTPUT ************************
1650	ccc************** WRF OUTPUT ************************
1651	ccc************** WRF OUTPUT ************************
1652	!do ig=1,ngrid
1653	! wtsurf(ig) = tsurf(ig) !! surface temperature
1654	! wco2ice(ig) = co2ice(ig) !! co2 ice
1655	!
1656	! !!! specific to WRF WRF WRF
1657	! !!! just to output water ice on surface
1658	! !!! uncomment the Registry entry
1659	! IF (igcm_h2o_ice .ne. 0) qsurfice(ig) = qsurf(ig,igcm_h2o_ice)
1660	!
1661	! !!! "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm"
1662	! IF (igcm_h2o_ice .ne. 0) THEN
1663	! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)*mugaz/mmol(igcm_h2o_ice)
1664	! ENDIF
1665	!
1666	!enddo
1667	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
1668	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
1669	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
1670	IF (igcm_h2o_ice .ne. 0) THEN
1671	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
1672	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1673	. * mugaz / mmol(igcm_h2o_ice)
1674	ENDIF
1675
1676	c
1677	c THIS INCLUDE IS AUTOMATICALLY GENERATED FROM REGISTRY
1678	c
1679	#include "fill_save.inc"
1680	c
1681	ccc************** WRF OUTPUT ************************
1682	ccc************** WRF OUTPUT ************************
1683	ccc************** WRF OUTPUT ************************
1684
1685
1686	cc-----------------------------------
1687	cc you can still use meso_WRITEDIAGFI (e.g. for debugging purpose),
1688	cc though this is not the default strategy now
1689	cc-----------------------------------
1690	cc please use cudt in namelist.input to set frequency of outputs
1691	cc-----------------------------------
1692	cc BEWARE: if at least one call to meso_WRITEDIAGFI is performed,
1693	cc cudt cannot be 0 - otherwise you'll get a "Floating exception"
1694	cc-----------------------------------
1695	! call meso_WRITEDIAGFI(ngrid,"tauref",
1696	! . "tauref","W.m-2",2,
1697	! . tauref)
1698	! call meso_WRITEDIAGFI(ngrid,"dtrad",
1699	! . "dtrad","W.m-2",2,
1700	! . dtrad)
1701	c call meso_WRITEDIAGFI(ngrid,"tsurf",
1702	c . "tsurf","K",2,
1703	c . tsurf)
1704	c
1705	! call meso_WRITEDIAGFI(ngrid,"zt",
1706	! . "zt","W.m-2",3,
1707	! . zt)
1708	! call meso_WRITEDIAGFI(ngrid,"zdtlw",
1709	! . "zdtlw","W.m-2",2,
1710	! . zdtlw)
1711	! call meso_WRITEDIAGFI(ngrid,"zdtsw",
1712	! . "zdtsw","W.m-2",2,
1713	! . zdtsw)
1714
1715
1716	!!
1717	!! ***WRF: everything below is kept for reference
1718	!!
1719	!
1720	!c ==========================================================
1721	!c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
1722	!c any variable for diagnostic (output with period
1723	!c "ecritphy", set in "run.def")
1724	!c ==========================================================
1725	!c WRITEDIAGFI can ALSO be called from any other subroutines
1726	!c for any variables !!
1727	! call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1728	! & emis)
1729	! call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1730	! & tsurf)
1731	! call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1732	! call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
1733	! & co2ice)
1734	!c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
1735	!c & "K",2,zt(1,7))
1736	! call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1737	! & fluxsurf_lw)
1738	! call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1739	! & fluxsurf_sw_tot)
1740	! call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1741	! & fluxtop_lw)
1742	! call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1743	! & fluxtop_sw_tot)
1744	! call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1745	!c call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1746	!c call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1747	!c call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1748	! call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1749	!c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1750	!c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1751	!c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
1752	!c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
1753	!c & zstress)
1754	!
1755	!c ----------------------------------------------------------
1756	!c Outputs of the CO2 cycle
1757	!c ----------------------------------------------------------
1758	!
1759	! if (tracer.and.(igcm_co2.ne.0)) then
1760	!! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
1761	!! & "kg/kg",2,zq(1,1,igcm_co2))
1762	!! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
1763	!! & "kg/kg",3,zq(1,1,igcm_co2))
1764	!
1765	! ! Compute co2 column
1766	! call zerophys(ngrid,co2col)
1767	! do l=1,nlayermx
1768	! do ig=1,ngrid
1769	! co2col(ig)=co2col(ig)+
1770	! & zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
1771	! enddo
1772	! enddo
1773	! call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
1774	! & co2col)
1775	! endif ! of if (tracer.and.(igcm_co2.ne.0))
1776	!
1777	!c ----------------------------------------------------------
1778	!c Outputs of the water cycle
1779	!c ----------------------------------------------------------
1780	! if (tracer) then
1781	! if (water) then
1782	!
1783	! CALL WRITEDIAGFI(ngridmx,'mtot',
1784	! & 'total mass of water vapor',
1785	! & 'kg/m2',2,mtot)
1786	! CALL WRITEDIAGFI(ngridmx,'icetot',
1787	! & 'total mass of water ice',
1788	! & 'kg/m2',2,icetot)
1789	!c If activice is true, tauTES is computed in aeropacity.F;
1790	! if (.not.activice) then
1791	! CALL WRITEDIAGFI(ngridmx,'tauTES',
1792	! & 'tau abs 825 cm-1',
1793	! & '',2,tauTES)
1794	! endif
1795	!
1796	! call WRITEDIAGFI(ngridmx,'h2o_ice_s',
1797	! & 'surface h2o_ice',
1798	! & 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
1799	!
1800	! if (activice) then
1801	!c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
1802	!c & 'w.m-2',3,zdtsw)
1803	!c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
1804	!c & 'w.m-2',3,zdtlw)
1805	! endif !(activice)
1806	! endif !(water)
1807	!
1808	!
1809	! if (water.and..not.photochem) then
1810	! iq=nq
1811	!c write(str2(1:2),'(i2.2)') iq
1812	!c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
1813	!c & 'kg.m-2',2,zdqscloud(1,iq))
1814	!c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
1815	!c & 'kg/kg',3,zdqchim(1,1,iq))
1816	!c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
1817	!c & 'kg/kg',3,zdqdif(1,1,iq))
1818	!c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
1819	!c & 'kg/kg',3,zdqadj(1,1,iq))
1820	!c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
1821	!c & 'kg/kg',3,zdqc(1,1,iq))
1822	! endif !(water.and..not.photochem)
1823	! endif
1824	!
1825	!c ----------------------------------------------------------
1826	!c Outputs of the dust cycle
1827	!c ----------------------------------------------------------
1828	!
1829	! call WRITEDIAGFI(ngridmx,'taudustvis',
1830	! & 'Dust optical depth',' ',2,tau(1,1))
1831	!
1832	! if (tracer.and.(dustbin.ne.0)) then
1833	! call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
1834	! if (doubleq) then
1835	! call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
1836	! & 'kg.m-2',2,qsurf(1,1))
1837	! call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
1838	! & 'N.m-2',2,qsurf(1,2))
1839	! call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
1840	! & 'kg.m-2.s-1',2,zdqsdev(1,1))
1841	! call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
1842	! & 'kg.m-2.s-1',2,zdqssed(1,1))
1843	! do l=1,nlayer
1844	! do ig=1, ngrid
1845	! reff(ig,l)= ref_r0 *
1846	! & (r3n_qpq(ig,l,1)/max(pq(ig,l,2),0.01))*(1./3.)
1847	! reff(ig,l)=min(max(reff(ig,l),1.e-10),500.e-6)
1848	! end do
1849	! end do
1850	! call WRITEDIAGFI(ngridmx,'reff','reff','m',3,reff)
1851	! else
1852	! do iq=1,dustbin
1853	! write(str2(1:2),'(i2.2)') iq
1854	! call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
1855	! & 'kg/kg',3,zq(1,1,iq))
1856	! call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
1857	! & 'kg.m-2',2,qsurf(1,iq))
1858	! end do
1859	! endif ! (doubleq)
1860	! end if ! (tracer.and.(dustbin.ne.0))
1861	!
1862	!c ----------------------------------------------------------
1863	!c Output in netcdf file "diagsoil.nc" for subterranean
1864	!c variables (output every "ecritphy", as for writediagfi)
1865	!c ----------------------------------------------------------
1866	!
1867	! ! Write soil temperature
1868	!! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
1869	!! & 3,tsoil)
1870	! ! Write surface temperature
1871	!! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
1872	!! & 2,tsurf)
1873	!
1874	!c ==========================================================
1875	!c END OF WRITEDIAGFI
1876	!c ==========================================================
1877
1878	ELSE ! if(ngrid.eq.1)
1879
1880	print,'Ls =',zls180./pi,
1881	& ' tauref(700 Pa) =',tauref
1882	c ----------------------------------------------------------------------
1883	c Output in grads file "g1d" (ONLY when using testphys1d)
1884	c (output at every X physical timestep)
1885	c ----------------------------------------------------------------------
1886	c
1887	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
1888	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
1889	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
1890
1891	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
1892	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
1893	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
1894	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
1895
1896	!! or output in diagfi.nc (for testphys1d)
1897	! call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
1898	! call WRITEDIAGFI(ngridmx,'temp','Temperature',
1899	! & 'K',1,zt)
1900	!
1901	! if(tracer) then
1902	!c CALL writeg1d(ngrid,1,tau,'tau','SI')
1903	! do iq=1,nq
1904	!c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
1905	! call WRITEDIAGFI(ngridmx,trim(noms(iq)),
1906	! & trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
1907	! end do
1908	! end if
1909	!
1910	! zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
1911	!
1912	! do l=2,nlayer-1
1913	! tmean=zt(1,l)
1914	! if(zt(1,l).ne.zt(1,l-1))
1915	! & tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
1916	! zlocal(l)= zlocal(l-1)
1917	! & -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
1918	! enddo
1919	! zlocal(nlayer)= zlocal(nlayer-1)-
1920	! & log(pplay(1,nlayer)/pplay(1,nlayer-1))*
1921	! & rnew(1,nlayer)*tmean/g
1922
1923	END IF ! if(ngrid.ne.1)
1924
1925	icount=icount+1
1926	write(,) 'now, back to the dynamical core...'
1927	#endif
1928	RETURN
1929	END

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

Download in other formats:

Original Format