Context Navigation

meso_physiq.F @ 40

Last change on this file since 40 was 35, checked in by aslmd, 15 years ago
LMD_MM_MARS: travail de preparation a un dossier libf commun avec LMDZ.MARS
File size: 64.0 KB

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

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

Download in other formats:

Original Format