Context Navigation

meso_physiq.F @ 161

Last change on this file since 161 was 131, checked in by aslmd, 14 years ago

mars

M 130 mars/libf/phymars/callradite.F
M 130 mars/README
no more need to modify callradite.F prior to compilation
[but still dimradmars.h must be modified]

--> in callradite.F we now have

-- DEFAULT name_iaer(1) is "dust_conrath"
-- IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq"
-- IF (water.AND.activice) name_iaer(2) = "h2o_ice"

M 130 000-USERS
small meeting about SVN for "planeto" team

M 130 mars/libf/phymars/meso_physiq.F
purely cosmetic changes

M 130 mesoscale/PLOT/MINIMAL/map_latlon.pro
M 130 mesoscale/PLOT/SPEC/GW/gravitwave2.pro
M 130 mesoscale/PLOT/SPEC/GW/gravitwaveprof.pro
M 130 mesoscale/PLOT/SPEC/GW/gravitwave_inc.pro
A 0 mesoscale/PLOT/SPEC/MAP/defs/THARSIS_newphys_bugnoRAC_RACgcm_TAUTES.map_uvt_inc.pro
M 130 mesoscale/PLOT/SPEC/MAP/map_uvt.pro
graphic stuff for GW studies and H2O cycle

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

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

Download in other formats:

Original Format