Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

physiq.F

Last change on this file was 1202, checked in by aslmd, 11 years ago
MESOSCALE. old physics for polar runs: better handling of Titus cap and handling of Tyler surface properties
File size: 51.2 KB

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

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

Download in other formats:

Original Format