Context Navigation

meso_physiq.F_old

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

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

Download in other formats:

Original Format