Context Navigation

meso_physiq.F @ 190

Last change on this file since 190 was 185, checked in by acolaitis, 13 years ago

17/06/2011 == AC

Added new settings for the Martian thermals from new LES observations
Revamped thermcell's module variables to allow it's removal
Minor changes in physiq and meso_physiq for the call to thermals
Switched from dynamic to static memory allocation for all thermals variable

to gain computation speed

Thermals now output maximum speed, maximum heat flux, and maximum height in thermal plumes

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

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

Download in other formats:

Original Format