Context Navigation

meso_physiq.F @ 61

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

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

Download in other formats:

Original Format