Context Navigation

meso_physiq.F @ 73

Last change on this file since 73 was 73, checked in by aslmd, 14 years ago
LMD_MM_MARS: corrections cycle de l'eau propagees a la nouvelle physique. + corrections readmeteo.F90 [version synchronisee precedemment n etait pas la plus a jour] + corrections api.F90 pour avoir cp, R comme GCM
File size: 65.3 KB

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

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

Download in other formats:

Original Format