Context Navigation

meso_physiq.F @ 91

Last change on this file since 91 was 86, checked in by aslmd, 14 years ago

*
mars + LMD_MM_MARS * Precompilation flag MESOSCALE for better transparency

* in shared phymars between GCM and mesoscale model

M 85 mars/libf/phymars/meso_physiq.F
M 85 mars/libf/phymars/meso_inifis.F
Added a pre-compilation flag MESOSCALE so that the LMDZ.MARS GCM
will compile without stating errors because of mesoscale routines.

M 85 mars/libf/phymars/newcondens.F
M 85 mars/libf/phymars/testphys1d.F
M 85 mars/libf/phymars/dustlift.F
D 85 mars/libf/phymars/meso_testphys1d.F
D 85 mars/libf/phymars/meso_dustlift.F
D 85 mars/libf/phymars/meso_newcondens.F
Now, this MESOSCALE precompilation flag can be used to lower
the number of meso_* routines when adaptations for mesoscale
applications are not very extended.
--> Three meso_* routines were deleted and changes are
now impacted under the MESOSCALE flag in the original GCM routines
--> Completely transparent for GCM compilation since it is devoid
of the -DMESOSCALE option
--> Very good for syncing because changes in dustlift, newcondens
will be directly available in the mesoscale model

M 84 mesoscale/LMD_MM_MARS/makemeso
Changed meso_testphys1d in testphys1d

M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_pgf
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpifort
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_ifort
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_g95
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpi
Added the option -DMESOSCALE in these scripts

*
LMD_MM_MARS * Various minor changes related to water cycle and plotting routines

* Also included the GW test case

A 0 mesoscale/LMDZ.MARS.new/myGCM/DEFS_JB/callphys.def.orig
M 84 mesoscale/NOTES.txt
D 84 mesoscale/LMD_MM_MARS/SRC/ARWpost/idl
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM
M 84 mesoscale/LMD_MM_MARS/SIMU/gnome_launch.meso
M 85 mesoscale/PLOT/MINIMAL/map_latlon.pro
D 85 mesoscale/PLOT/SPEC/LES/getget.pro
M 85 mesoscale/PLOT/SPEC/MAP/map_uvt.pro
A + - mesoscale/PLOT/SPEC/getget.pro
A 0 mesoscale/PLOT/RESERVE/obsolete
A 0 mesoscale/TESTS/TESTGW.tar.gz
M 84 000-USERS

File size: 67.8 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	!!!**WRF: newphys: ici la pression n'est plus mise a jour ds le GCM
1063	!!!**WRF: mais il faut retablir ca dans le cas du mesoscale ?
1064	!!!**WRF: ...non probablement OK
1065	! ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1066	ENDDO
1067
1068	IF (tracer) THEN
1069	DO iq=1, nq
1070	DO l=1,nlayer
1071	DO ig=1,ngrid
1072	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1073	ENDDO
1074	ENDDO
1075	ENDDO
1076	ENDIF ! of IF (tracer)
1077
1078	ENDIF ! of IF (callcond)
1079
1080	c-----------------------------------------------------------------------
1081	c 7. Specific parameterizations for tracers
1082	c: -----------------------------------------
1083
1084	if (tracer) then
1085
1086	c 7a. Water and ice
1087	c ---------------
1088
1089	c ---------------------------------------
1090	c Water ice condensation in the atmosphere
1091	c ----------------------------------------
1092	IF (water) THEN
1093
1094	c **WRF: new arguments here rnuclei,rice,nuice
1095	c plus no more iqmin +igcm_h2o_vap replaces iq, what are the consequences?
1096	c checks needed when tracers simulations
1097
1098	call watercloud(ngrid,nlayer,ptimestep,
1099	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
1100	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
1101	& nq,naerkind,tau,
1102	& ccn,rdust,rice,nuice)
1103	if (activice) then
1104	c Temperature variation due to latent heat release
1105	DO l=1,nlayer
1106	DO ig=1,ngrid
1107	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
1108	ENDDO
1109	ENDDO
1110	endif
1111
1112	! increment water vapour and ice atmospheric tracers tendencies
1113	IF (water) THEN
1114	DO l=1,nlayer
1115	DO ig=1,ngrid
1116	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
1117	& zdqcloud(ig,l,igcm_h2o_vap)
1118	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
1119	& zdqcloud(ig,l,igcm_h2o_ice)
1120	ENDDO
1121	ENDDO
1122	ENDIF ! of IF (water) THEN
1123	! Increment water ice surface tracer tendency
1124	DO ig=1,ngrid
1125	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
1126	& zdqscloud(ig,igcm_h2o_ice)
1127	ENDDO
1128
1129	END IF ! of IF (water)
1130
1131	c 7b. Chemical species
1132	c ------------------
1133
1134	!!!
1135	!!! WRF WRF WRF commented for smaller executables
1136	!!!
1137	!c --------------
1138	!c photochemistry :
1139	!c --------------
1140	! IF (photochem .or. thermochem) then
1141	! call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
1142	! & zzlay,zday,pq,pdq,rice,
1143	! & zdqchim,zdqschim,zdqcloud,zdqscloud)
1144	!!NB: Photochemistry includes condensation of H2O2
1145	!
1146	! ! increment values of tracers:
1147	! DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1148	! ! tracers is zero anyways
1149	! DO l=1,nlayer
1150	! DO ig=1,ngrid
1151	! pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
1152	! ENDDO
1153	! ENDDO
1154	! ENDDO ! of DO iq=1,nq
1155	! ! add condensation tendency for H2O2
1156	! if (igcm_h2o2.ne.0) then
1157	! DO l=1,nlayer
1158	! DO ig=1,ngrid
1159	! pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
1160	! & +zdqcloud(ig,l,igcm_h2o2)
1161	! ENDDO
1162	! ENDDO
1163	! endif
1164	!
1165	! ! increment surface values of tracers:
1166	! DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1167	! ! tracers is zero anyways
1168	! DO ig=1,ngrid
1169	! dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
1170	! ENDDO
1171	! ENDDO ! of DO iq=1,nq
1172	! ! add condensation tendency for H2O2
1173	! if (igcm_h2o2.ne.0) then
1174	! DO ig=1,ngrid
1175	! dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
1176	! & +zdqscloud(ig,igcm_h2o2)
1177	! ENDDO
1178	! endif
1179	!
1180	! END IF ! of IF (photochem.or.thermochem)
1181
1182	c 7c. Aerosol particles
1183	c -------------------
1184
1185	c ----------
1186	c Dust devil :
1187	c ----------
1188	IF(callddevil) then
1189	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
1190	& zdqdev,zdqsdev)
1191
1192	if (dustbin.ge.1) then
1193	do iq=1,nq
1194	DO l=1,nlayer
1195	DO ig=1,ngrid
1196	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
1197	ENDDO
1198	ENDDO
1199	enddo
1200	do iq=1,nq
1201	DO ig=1,ngrid
1202	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
1203	ENDDO
1204	enddo
1205	endif ! of if (dustbin.ge.1)
1206
1207	END IF ! of IF (callddevil)
1208
1209	c -------------
1210	c Sedimentation : acts also on water ice
1211	c -------------
1212	IF (sedimentation) THEN
1213	!call zerophys(ngridnlayernq, zdqsed)
1214	zdqsed(1:ngrid,1:nlayer,1:nq)=0
1215	!call zerophys(ngrid*nq, zdqssed)
1216	zdqssed(1:ngrid,1:nq)=0
1217
1218	c
1219	c **WRF: new arguments rnuclei, rice, need checks
1220	c
1221	call callsedim(ngrid,nlayer, ptimestep,
1222	& pplev,zzlev, pt, rdust, rice,
1223	& pq, pdq, zdqsed, zdqssed,nq)
1224
1225	DO iq=1, nq
1226	DO l=1,nlayer
1227	DO ig=1,ngrid
1228	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1229	ENDDO
1230	ENDDO
1231	ENDDO
1232	DO iq=1, nq
1233	DO ig=1,ngrid
1234	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1235	ENDDO
1236	ENDDO
1237	END IF ! of IF (sedimentation)
1238
1239	c 7d. Updates
1240	c ---------
1241
1242	DO iq=1, nq
1243	DO ig=1,ngrid
1244
1245	c ---------------------------------
1246	c Updating tracer budget on surface
1247	c ---------------------------------
1248	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1249
1250	ENDDO ! (ig)
1251	ENDDO ! (iq)
1252
1253	endif ! of if (tracer)
1254
1255	!!!
1256	!!! WRF WRF WRF commented for smaller executables
1257	!!!
1258	!c-----------------------------------------------------------------------
1259	!c 8. THERMOSPHERE CALCULATION
1260	!c-----------------------------------------------------------------------
1261	!
1262	! if (callthermos) then
1263	! call thermosphere(pplev,pplay,dist_sol,
1264	! $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1265	! & pt,pq,pu,pv,pdt,pdq,
1266	! $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1267	!
1268	! DO l=1,nlayer
1269	! DO ig=1,ngrid
1270	! dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1271	! pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1272	! & +zdteuv(ig,l)
1273	! pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1274	! pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1275	! DO iq=1, nq
1276	! pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1277	! ENDDO
1278	! ENDDO
1279	! ENDDO
1280	!
1281	! endif ! of if (callthermos)
1282
1283	c-----------------------------------------------------------------------
1284	c 9. Surface and sub-surface soil temperature
1285	c-----------------------------------------------------------------------
1286	c
1287	c
1288	c 9.1 Increment Surface temperature:
1289	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1290
1291	DO ig=1,ngrid
1292	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1293	ENDDO
1294
1295	ccc
1296	ccc **WRF very specific to GCM
1297	ccc
1298	c Prescribe a cold trap at south pole (except at high obliquity !!)
1299	c Temperature at the surface is set there to be the temperature
1300	c corresponding to equilibrium temperature between phases of CO2
1301
1302	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
1303	! if (caps.and.(obliquit.lt.27.)) then
1304	! ! NB: Updated surface pressure, at grid point 'ngrid', is
1305	! ! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1306	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
1307	! & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1308	! endif
1309	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1310	!!!!! note WRF MESOSCALE AYMERIC -- mot cle "caps"
1311	!!!!! watercaptag n'est plus utilise que dans vdifc
1312	!!!!! ... pour que la sublimation ne soit pas stoppee
1313	!!!!! ... dans la calotte permanente nord si qsurf=0
1314	!!!!! on desire garder cet effet regle par caps=T
1315	!!!!! on a donc commente "if (caps.and.(obliquit.lt.27.))" ci-dessus
1316	!!!!! --- remplacer ces lignes par qqch de plus approprie
1317	!!!!! si on s attaque a la calotte polaire sud
1318	!!!!! pas d'autre occurrence majeure du mot-cle "caps"
1319	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1320
1321	c -------------------------------------------------------------
1322	c Change of surface albedo (set to 0.4) in case of ground frost
1323	c everywhere except on the north permanent cap and in regions
1324	c covered by dry ice.
1325	c ALWAYS PLACE these lines after newcondens !!!
1326	c -------------------------------------------------------------
1327	c
1328	c **WRF : OK avec le mesoscale, pas d'indices bizarres au pole
1329	c
1330	do ig=1,ngrid
1331	if ((co2ice(ig).eq.0).and.
1332	& (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
1333	albedo(ig,1) = alb_surfice
1334	albedo(ig,2) = alb_surfice
1335	endif
1336	enddo ! of do ig=1,ngrid
1337	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
1338
1339	c
1340	c 9.2 Compute soil temperatures and subsurface heat flux:
1341	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1342	IF (callsoil) THEN
1343	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1344	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
1345	ENDIF
1346
1347	c-----------------------------------------------------------------------
1348	c 10. Write output files
1349	c ----------------------
1350
1351	c -------------------------------
1352	c Dynamical fields incrementation
1353	c -------------------------------
1354	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1355	! temperature, zonal and meridional wind
1356	DO l=1,nlayer
1357	DO ig=1,ngrid
1358	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
1359	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
1360	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
1361	ENDDO
1362	ENDDO
1363
1364	! tracers
1365	DO iq=1, nq
1366	DO l=1,nlayer
1367	DO ig=1,ngrid
1368	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1369	ENDDO
1370	ENDDO
1371	ENDDO
1372
1373	! surface pressure
1374	DO ig=1,ngrid
1375	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
1376	ENDDO
1377
1378	! pressure
1379	DO l=1,nlayer
1380	DO ig=1,ngrid
1381	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1382	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1383	ENDDO
1384	ENDDO
1385
1386	! Density
1387	DO l=1,nlayer
1388	DO ig=1,ngrid
1389	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1390	ENDDO
1391	ENDDO
1392
1393	c Compute surface stress : (NB: z0 is a common in planete.h)
1394	c DO ig=1,ngrid
1395	c cd = (0.4/log(zzlay(ig,1)/z0))**2
1396	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
1397	c ENDDO
1398
1399	c Sum of fluxes in solar spectral bands (for output only)
1400	DO ig=1,ngrid
1401	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1402	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1403	ENDDO
1404	c ***** TEST ****************************************************
1405	ztim1 = 999
1406	DO l=1,nlayer
1407	DO ig=1,ngrid
1408	if (pt(ig,l).lt.ztim1) then
1409	ztim1 = pt(ig,l)
1410	igmin = ig
1411	lmin = l
1412	end if
1413	ENDDO
1414	ENDDO
1415	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1416	write(,) 'PHYSIQ: stability WARNING :'
1417	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1418	& 'ig l =', igmin, lmin
1419	end if
1420	c *******************************************************************
1421
1422	c ---------------------
1423	c Outputs to the screen
1424	c ---------------------
1425
1426	IF (lwrite) THEN
1427	PRINT*,'Global diagnostics for the physics'
1428	PRINT,'Variables and their increments x and dx/dt dt'
1429	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1430	WRITE(*,'(2f10.5,2f15.5)')
1431	s tsurf(igout),zdtsurf(igout)*ptimestep,
1432	s pplev(igout,1),pdpsrf(igout)*ptimestep
1433	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1434	WRITE(*,'(i4,6f10.5)') (l,
1435	s pu(igout,l),pdu(igout,l)*ptimestep,
1436	s pv(igout,l),pdv(igout,l)*ptimestep,
1437	s pt(igout,l),pdt(igout,l)*ptimestep,
1438	s l=1,nlayer)
1439	ENDIF ! of IF (lwrite)
1440
1441	IF (ngrid.NE.1) THEN
1442	print,'Ls =',zls180./pi,
1443	& ' tauref(700 Pa,lat=0) =',tauref(ngrid/2)!,
1444	! & ' tau(Viking1) =',tau(ig_vl1,1)
1445
1446
1447	c -------------------------------------------------------------------
1448	c Writing NetCDF file "RESTARTFI" at the end of the run
1449	c -------------------------------------------------------------------
1450	c Note: 'restartfi' is stored just before dynamics are stored
1451	c in 'restart'. Between now and the writting of 'restart',
1452	c there will have been the itau=itau+1 instruction and
1453	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1454	c thus we store for time=time+dtvr
1455
1456
1457	!!!
1458	!!! WRF WRF WRF WRF
1459	!!!
1460	! IF(lastcall) THEN
1461	! ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1462	! write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1463	! call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
1464	! . ptimestep,pday,
1465	! . ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
1466	! . area,albedodat,inertiedat,zmea,zstd,zsig,
1467	! . zgam,zthe)
1468	! ENDIF
1469
1470
1471
1472	c -------------------------------------------------------------------
1473	c Calculation of diagnostic variables written in both stats and
1474	c diagfi files
1475	c -------------------------------------------------------------------
1476
1477	if (tracer) then
1478	if (water) then
1479
1480	!!
1481	!!***WRF: ok, des nouveaux trucs cools de la nouvelle physique
1482	!!
1483	call zerophys(ngrid,mtot)
1484	call zerophys(ngrid,icetot)
1485	call zerophys(ngrid,rave)
1486	call zerophys(ngrid,tauTES)
1487	do ig=1,ngrid
1488	do l=1,nlayermx
1489	mtot(ig) = mtot(ig) +
1490	& zq(ig,l,igcm_h2o_vap) *
1491	& (pplev(ig,l) - pplev(ig,l+1)) / g
1492	icetot(ig) = icetot(ig) +
1493	& zq(ig,l,igcm_h2o_ice) *
1494	& (pplev(ig,l) - pplev(ig,l+1)) / g
1495	rave(ig) = rave(ig) +
1496	& zq(ig,l,igcm_h2o_ice) *
1497	& (pplev(ig,l) - pplev(ig,l+1)) / g *
1498	& rice(ig,l) * (1.+nuice_ref)
1499	c Computing abs optical depth at 825 cm-1 in each
1500	c layer to simulate NEW TES retrieval
1501	Qabsice = min(
1502	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
1503	& )
1504	opTES(ig,l)= 0.75 * Qabsice *
1505	& zq(ig,l,igcm_h2o_ice) *
1506	& (pplev(ig,l) - pplev(ig,l+1)) / g
1507	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1508	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1509	enddo
1510	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
1511	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1512	enddo
1513
1514	endif ! of if (water)
1515	endif ! of if (tracer)
1516
1517	c -----------------------------------------------------------------
1518	c WSTATS: Saving statistics
1519	c -----------------------------------------------------------------
1520	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
1521	c which can later be used to make the statistic files of the run:
1522	c "stats") only possible in 3D runs !
1523
1524
1525	IF (callstats) THEN
1526
1527	write(,) 'callstats'
1528
1529	! call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1530	! call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1531	! call wstats(ngrid,"co2ice","CO2 ice cover",
1532	! & "kg.m-2",2,co2ice)
1533	! call wstats(ngrid,"fluxsurf_lw",
1534	! & "Thermal IR radiative flux to surface","W.m-2",2,
1535	! & fluxsurf_lw)
1536	! call wstats(ngrid,"fluxsurf_sw",
1537	! & "Solar radiative flux to surface","W.m-2",2,
1538	! & fluxsurf_sw_tot)
1539	! call wstats(ngrid,"fluxtop_lw",
1540	! & "Thermal IR radiative flux to space","W.m-2",2,
1541	! & fluxtop_lw)
1542	! call wstats(ngrid,"fluxtop_sw",
1543	! & "Solar radiative flux to space","W.m-2",2,
1544	! & fluxtop_sw_tot)
1545	! call wstats(ngrid,"taudustvis",
1546	! & "Dust optical depth"," ",2,tau(1,1))
1547	! call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1548	! call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1549	! call wstats(ngrid,"v","Meridional (North-South) wind",
1550	! & "m.s-1",3,zv)
1551	!c call wstats(ngrid,"w","Vertical (down-up) wind",
1552	!c & "m.s-1",3,pw)
1553	! call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
1554	!c call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
1555	!c call wstats(ngrid,"q2",
1556	!c & "Boundary layer eddy kinetic energy",
1557	!c & "m2.s-2",3,q2)
1558	!c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1559	!c & emis)
1560	!c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
1561	!c & 2,zstress)
1562	!
1563	! if (tracer) then
1564	! if (water) then
1565	! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
1566	! & *mugaz/mmol(igcm_h2o_vap)
1567	! call wstats(ngrid,"vmr_h2ovapor",
1568	! & "H2O vapor volume mixing ratio","mol/mol",
1569	! & 3,vmr)
1570	! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
1571	! & *mugaz/mmol(igcm_h2o_ice)
1572	! call wstats(ngrid,"vmr_h2oice",
1573	! & "H2O ice volume mixing ratio","mol/mol",
1574	! & 3,vmr)
1575	!
1576	! call wstats(ngrid,"mtot",
1577	! & "total mass of water vapor","kg/m2",
1578	! & 2,mtot)
1579	! call wstats(ngrid,"icetot",
1580	! & "total mass of water ice","kg/m2",
1581	! & 2,icetot)
1582	!c If activice is true, tauTES is computed in aeropacity.F;
1583	! if (.not.activice) then
1584	! call wstats(ngrid,"tauTES",
1585	! & "tau abs 825 cm-1","",
1586	! & 2,tauTES)
1587	! endif ! of if (activice)
1588	!
1589	! endif ! of if (water)
1590	!
1591	! if (thermochem.or.photochem) then
1592	! do iq=1,nq
1593	! if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
1594	! . (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
1595	! . (noms(iq).eq."h2").or.
1596	! . (noms(iq).eq."o3")) then
1597	! do l=1,nlayer
1598	! do ig=1,ngrid
1599	! vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
1600	! end do
1601	! end do
1602	! call wstats(ngrid,"vmr_"//trim(noms(iq)),
1603	! . "Volume mixing ratio","mol/mol",3,vmr)
1604	! endif
1605	! enddo
1606	! endif ! of if (thermochem.or.photochem)
1607	!
1608	! endif ! of if (tracer)
1609	!
1610	! IF(lastcall) THEN
1611	! write (,) "Writing stats..."
1612	! call mkstats(ierr)
1613	! ENDIF
1614
1615	ENDIF !if callstats
1616
1617	c (Store EOF for Mars Climate database software)
1618	IF (calleofdump) THEN
1619	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1620	ENDIF
1621
1622	ccc************** WRF OUTPUT ************************
1623	ccc************** WRF OUTPUT ************************
1624	ccc************** WRF OUTPUT ************************
1625	!do ig=1,ngrid
1626	! wtsurf(ig) = tsurf(ig) !! surface temperature
1627	! wco2ice(ig) = co2ice(ig) !! co2 ice
1628	!
1629	! !!! specific to WRF WRF WRF
1630	! !!! just to output water ice on surface
1631	! !!! uncomment the Registry entry
1632	! IF (igcm_h2o_ice .ne. 0) qsurfice(ig) = qsurf(ig,igcm_h2o_ice)
1633	!
1634	! !!! "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm"
1635	! IF (igcm_h2o_ice .ne. 0) THEN
1636	! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)*mugaz/mmol(igcm_h2o_ice)
1637	! ENDIF
1638	!
1639	!enddo
1640	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
1641	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
1642	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
1643	IF (igcm_h2o_ice .ne. 0) THEN
1644	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
1645	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1646	. * mugaz / mmol(igcm_h2o_ice)
1647	ENDIF
1648
1649	c
1650	c THIS INCLUDE IS AUTOMATICALLY GENERATED FROM REGISTRY
1651	c
1652	#include "fill_save.inc"
1653	c
1654	ccc************** WRF OUTPUT ************************
1655	ccc************** WRF OUTPUT ************************
1656	ccc************** WRF OUTPUT ************************
1657
1658
1659	cc-----------------------------------
1660	cc you can still use meso_WRITEDIAGFI (e.g. for debugging purpose),
1661	cc though this is not the default strategy now
1662	cc-----------------------------------
1663	cc please use cudt in namelist.input to set frequency of outputs
1664	cc-----------------------------------
1665	cc BEWARE: if at least one call to meso_WRITEDIAGFI is performed,
1666	cc cudt cannot be 0 - otherwise you'll get a "Floating exception"
1667	cc-----------------------------------
1668	! call meso_WRITEDIAGFI(ngrid,"tauref",
1669	! . "tauref","W.m-2",2,
1670	! . tauref)
1671	! call meso_WRITEDIAGFI(ngrid,"dtrad",
1672	! . "dtrad","W.m-2",2,
1673	! . dtrad)
1674	c call meso_WRITEDIAGFI(ngrid,"tsurf",
1675	c . "tsurf","K",2,
1676	c . tsurf)
1677	c
1678	! call meso_WRITEDIAGFI(ngrid,"zt",
1679	! . "zt","W.m-2",3,
1680	! . zt)
1681	! call meso_WRITEDIAGFI(ngrid,"zdtlw",
1682	! . "zdtlw","W.m-2",2,
1683	! . zdtlw)
1684	! call meso_WRITEDIAGFI(ngrid,"zdtsw",
1685	! . "zdtsw","W.m-2",2,
1686	! . zdtsw)
1687
1688
1689	!!
1690	!! ***WRF: everything below is kept for reference
1691	!!
1692	!
1693	!c ==========================================================
1694	!c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
1695	!c any variable for diagnostic (output with period
1696	!c "ecritphy", set in "run.def")
1697	!c ==========================================================
1698	!c WRITEDIAGFI can ALSO be called from any other subroutines
1699	!c for any variables !!
1700	! call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1701	! & emis)
1702	! call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1703	! & tsurf)
1704	! call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1705	! call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
1706	! & co2ice)
1707	!c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
1708	!c & "K",2,zt(1,7))
1709	! call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1710	! & fluxsurf_lw)
1711	! call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1712	! & fluxsurf_sw_tot)
1713	! call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1714	! & fluxtop_lw)
1715	! call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1716	! & fluxtop_sw_tot)
1717	! call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1718	!c call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1719	!c call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1720	!c call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1721	! call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1722	!c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1723	!c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1724	!c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
1725	!c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
1726	!c & zstress)
1727	!
1728	!c ----------------------------------------------------------
1729	!c Outputs of the CO2 cycle
1730	!c ----------------------------------------------------------
1731	!
1732	! if (tracer.and.(igcm_co2.ne.0)) then
1733	!! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
1734	!! & "kg/kg",2,zq(1,1,igcm_co2))
1735	!! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
1736	!! & "kg/kg",3,zq(1,1,igcm_co2))
1737	!
1738	! ! Compute co2 column
1739	! call zerophys(ngrid,co2col)
1740	! do l=1,nlayermx
1741	! do ig=1,ngrid
1742	! co2col(ig)=co2col(ig)+
1743	! & zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
1744	! enddo
1745	! enddo
1746	! call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
1747	! & co2col)
1748	! endif ! of if (tracer.and.(igcm_co2.ne.0))
1749	!
1750	!c ----------------------------------------------------------
1751	!c Outputs of the water cycle
1752	!c ----------------------------------------------------------
1753	! if (tracer) then
1754	! if (water) then
1755	!
1756	! CALL WRITEDIAGFI(ngridmx,'mtot',
1757	! & 'total mass of water vapor',
1758	! & 'kg/m2',2,mtot)
1759	! CALL WRITEDIAGFI(ngridmx,'icetot',
1760	! & 'total mass of water ice',
1761	! & 'kg/m2',2,icetot)
1762	!c If activice is true, tauTES is computed in aeropacity.F;
1763	! if (.not.activice) then
1764	! CALL WRITEDIAGFI(ngridmx,'tauTES',
1765	! & 'tau abs 825 cm-1',
1766	! & '',2,tauTES)
1767	! endif
1768	!
1769	! call WRITEDIAGFI(ngridmx,'h2o_ice_s',
1770	! & 'surface h2o_ice',
1771	! & 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
1772	!
1773	! if (activice) then
1774	!c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
1775	!c & 'w.m-2',3,zdtsw)
1776	!c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
1777	!c & 'w.m-2',3,zdtlw)
1778	! endif !(activice)
1779	! endif !(water)
1780	!
1781	!
1782	! if (water.and..not.photochem) then
1783	! iq=nq
1784	!c write(str2(1:2),'(i2.2)') iq
1785	!c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
1786	!c & 'kg.m-2',2,zdqscloud(1,iq))
1787	!c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
1788	!c & 'kg/kg',3,zdqchim(1,1,iq))
1789	!c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
1790	!c & 'kg/kg',3,zdqdif(1,1,iq))
1791	!c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
1792	!c & 'kg/kg',3,zdqadj(1,1,iq))
1793	!c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
1794	!c & 'kg/kg',3,zdqc(1,1,iq))
1795	! endif !(water.and..not.photochem)
1796	! endif
1797	!
1798	!c ----------------------------------------------------------
1799	!c Outputs of the dust cycle
1800	!c ----------------------------------------------------------
1801	!
1802	! call WRITEDIAGFI(ngridmx,'taudustvis',
1803	! & 'Dust optical depth',' ',2,tau(1,1))
1804	!
1805	! if (tracer.and.(dustbin.ne.0)) then
1806	! call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
1807	! if (doubleq) then
1808	! call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
1809	! & 'kg.m-2',2,qsurf(1,1))
1810	! call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
1811	! & 'N.m-2',2,qsurf(1,2))
1812	! call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
1813	! & 'kg.m-2.s-1',2,zdqsdev(1,1))
1814	! call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
1815	! & 'kg.m-2.s-1',2,zdqssed(1,1))
1816	! do l=1,nlayer
1817	! do ig=1, ngrid
1818	! reff(ig,l)= ref_r0 *
1819	! & (r3n_qpq(ig,l,1)/max(pq(ig,l,2),0.01))*(1./3.)
1820	! reff(ig,l)=min(max(reff(ig,l),1.e-10),500.e-6)
1821	! end do
1822	! end do
1823	! call WRITEDIAGFI(ngridmx,'reff','reff','m',3,reff)
1824	! else
1825	! do iq=1,dustbin
1826	! write(str2(1:2),'(i2.2)') iq
1827	! call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
1828	! & 'kg/kg',3,zq(1,1,iq))
1829	! call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
1830	! & 'kg.m-2',2,qsurf(1,iq))
1831	! end do
1832	! endif ! (doubleq)
1833	! end if ! (tracer.and.(dustbin.ne.0))
1834	!
1835	!c ----------------------------------------------------------
1836	!c Output in netcdf file "diagsoil.nc" for subterranean
1837	!c variables (output every "ecritphy", as for writediagfi)
1838	!c ----------------------------------------------------------
1839	!
1840	! ! Write soil temperature
1841	!! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
1842	!! & 3,tsoil)
1843	! ! Write surface temperature
1844	!! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
1845	!! & 2,tsurf)
1846	!
1847	!c ==========================================================
1848	!c END OF WRITEDIAGFI
1849	!c ==========================================================
1850
1851	ELSE ! if(ngrid.eq.1)
1852
1853	print,'Ls =',zls180./pi,
1854	& ' tauref(700 Pa) =',tauref
1855	c ----------------------------------------------------------------------
1856	c Output in grads file "g1d" (ONLY when using testphys1d)
1857	c (output at every X physical timestep)
1858	c ----------------------------------------------------------------------
1859	c
1860	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
1861	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
1862	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
1863
1864	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
1865	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
1866	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
1867	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
1868
1869	!! or output in diagfi.nc (for testphys1d)
1870	! call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
1871	! call WRITEDIAGFI(ngridmx,'temp','Temperature',
1872	! & 'K',1,zt)
1873	!
1874	! if(tracer) then
1875	!c CALL writeg1d(ngrid,1,tau,'tau','SI')
1876	! do iq=1,nq
1877	!c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
1878	! call WRITEDIAGFI(ngridmx,trim(noms(iq)),
1879	! & trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
1880	! end do
1881	! end if
1882	!
1883	! zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
1884	!
1885	! do l=2,nlayer-1
1886	! tmean=zt(1,l)
1887	! if(zt(1,l).ne.zt(1,l-1))
1888	! & tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
1889	! zlocal(l)= zlocal(l-1)
1890	! & -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
1891	! enddo
1892	! zlocal(nlayer)= zlocal(nlayer-1)-
1893	! & log(pplay(1,nlayer)/pplay(1,nlayer-1))*
1894	! & rnew(1,nlayer)*tmean/g
1895
1896	END IF ! if(ngrid.ne.1)
1897
1898	icount=icount+1
1899	write(,) 'now, back to the dynamical core...'
1900	#endif
1901	RETURN
1902	END

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

Download in other formats:

Original Format