Context Navigation

meso_physiq.F @ 54

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

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

Download in other formats:

Original Format