Context Navigation

vdifc_mod.F @ 2953

Last change on this file since 2953 was 2953, checked in by romain.vande, 19 months ago
Mars PCM : Adapt vdifc, co2condens, rocketduststorm and topmons routines to the subslope parametrisation. RV
File size: 43.7 KB

Line
1	MODULE vdifc_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	SUBROUTINE vdifc(ngrid,nlay,nq,ppopsk,
8	$ ptimestep,pcapcal,lecrit,
9	$ pplay,pplev,pzlay,pzlev,pz0,
10	$ pu,pv,ph,pq,ptsrf,pemis,pqsurf,
11	$ pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
12	$ pdudif,pdvdif,pdhdif,pdtsrf,pq2,
13	$ pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
14	$ hfmax,pcondicea_co2microp,sensibFlux,
15	$ dustliftday,local_time,watercap, dwatercap_dif)
16
17	use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
18	& igcm_dust_submicron, igcm_h2o_vap,
19	& igcm_h2o_ice, alpha_lift, igcm_co2,
20	& igcm_hdo_vap, igcm_hdo_ice,
21	& igcm_stormdust_mass, igcm_stormdust_number
22	use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness
23	USE comcstfi_h, ONLY: cpp, r, rcp, g, pi
24	use watersat_mod, only: watersat
25	use turb_mod, only: turb_resolved, ustar, tstar
26	use compute_dtau_mod, only: ti_injection_sol,tf_injection_sol
27	use hdo_surfex_mod, only: hdo_surfex
28	c use geometry_mod, only: longitude_deg,latitude_deg ! Joseph
29	use dust_param_mod, only: doubleq, submicron, lifting
30	use write_output_mod, only: write_output
31	use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
32	& subslope_dist,major_slope,iflat
33
34	IMPLICIT NONE
35
36	c=======================================================================
37	c
38	c subject:
39	c --------
40	c Turbulent diffusion (mixing) for potential T, U, V and tracer
41	c
42	c Shema implicite
43	c On commence par rajouter au variables x la tendance physique
44	c et on resoult en fait:
45	c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
46	c
47	c author:
48	c ------
49	c Hourdin/Forget/Fournier
50	c=======================================================================
51
52	c-----------------------------------------------------------------------
53	c declarations:
54	c -------------
55
56	include "callkeys.h"
57	include "microphys.h"
58
59	c
60	c arguments:
61	c ----------
62
63	INTEGER,INTENT(IN) :: ngrid,nlay
64	REAL,INTENT(IN) :: ptimestep
65	REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
66	REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
67	REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
68	REAL,INTENT(IN) :: ph(ngrid,nlay)
69	REAL,INTENT(IN) :: ptsrf(ngrid,nslope),pemis(ngrid,nslope)
70	REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
71	REAL,INTENT(IN) :: pdhfi(ngrid,nlay)
72	REAL,INTENT(IN) :: pfluxsrf(ngrid,nslope)
73	REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
74	REAL,INTENT(OUT) :: pdtsrf(ngrid,nslope),pdhdif(ngrid,nlay)
75	REAL,INTENT(IN) :: pcapcal(ngrid,nslope)
76	REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
77
78	c Argument added for condensation:
79	REAL,INTENT(IN) :: ppopsk(ngrid,nlay)
80	logical,INTENT(IN) :: lecrit
81	REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlay)! tendency due to CO2 condensation (kg/kg.s-1)
82
83	REAL,INTENT(IN) :: pz0(ngrid) ! surface roughness length (m)
84
85	c Argument added to account for subgrid gustiness :
86
87	REAL,INTENT(IN) :: wstar(ngrid), hfmax(ngrid)!, zi(ngrid)
88
89	c Traceurs :
90	integer,intent(in) :: nq
91	REAL,INTENT(IN) :: pqsurf(ngrid,nq,nslope)
92	REAL :: zqsurf(ngrid) ! temporary water tracer
93	real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
94	real,intent(out) :: pdqdif(ngrid,nlay,nq)
95	real,intent(out) :: pdqsdif(ngrid,nq,nslope)
96	REAL,INTENT(in) :: dustliftday(ngrid)
97	REAL,INTENT(in) :: local_time(ngrid)
98
99	c local:
100	c ------
101
102	REAL :: pt(ngrid,nlay)
103
104	INTEGER ilev,ig,ilay,nlev,islope
105
106	REAL z4st,zdplanck(ngrid)
107	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
108	REAL zkq(ngrid,nlay+1)
109	REAL zcdv(ngrid),zcdh(ngrid)
110	REAL, INTENT(IN) :: zcdv_true(ngrid),zcdh_true(ngrid) ! drag coeff are used by the LES to recompute u* and hfx
111	REAL zu(ngrid,nlay),zv(ngrid,nlay)
112	REAL zh(ngrid,nlay)
113	REAL ztsrf2(ngrid)
114	REAL z1(ngrid),z2(ngrid)
115	REAL za(ngrid,nlay),zb(ngrid,nlay)
116	REAL zb0(ngrid,nlay)
117	REAL zc(ngrid,nlay),zd(ngrid,nlay)
118	REAL zcst1
119	REAL zu2(ngrid)
120
121	EXTERNAL SSUM,SCOPY
122	REAL SSUM
123	LOGICAL,SAVE :: firstcall=.true.
124
125	!$OMP THREADPRIVATE(firstcall)
126
127	c variable added for CO2 condensation:
128	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
129	REAL hh , zhcond(ngrid,nlay)
130	REAL,PARAMETER :: latcond=5.9e5
131	REAL,PARAMETER :: tcond1mb=136.27
132	REAL,SAVE :: acond,bcond
133
134	!$OMP THREADPRIVATE(acond,bcond)
135
136	c Subtimestep & implicit treatment of water vapor
137	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
138	REAL zdqsdif(ngrid) ! subtimestep pdqsdif for water ice
139	REAL ztsrf(ngrid) ! temporary surface temperature in tsub
140	REAL zdtsrf(ngrid,nslope) ! surface temperature tendancy in tsub
141	REAL surf_h2o_lh(ngrid,nslope) ! Surface h2o latent heat flux
142	REAL zsurf_h2o_lh(ngrid,nslope) ! Tsub surface h2o latent heat flux
143
144	c For latent heat release from ground water ice sublimation
145	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
146	REAL tsrf_lh(ngrid) ! temporary surface temperature with lh effect
147	REAL lh ! latent heat, formulation given in the Technical Document:
148	! "Modeling water ice sublimation under Phoenix-like conditions", Montmessin et al. 2004
149
150	c Tracers :
151	c ~~~~~~~
152	INTEGER iq
153	REAL zq(ngrid,nlay,nq)
154	REAL zq1temp(ngrid)
155	REAL rho(ngrid) ! near surface air density
156	REAL qsat(ngrid)
157
158	REAL hdoflux(ngrid,nslope) ! value of vapour flux of HDO
159	REAL hdoflux_meshavg(ngrid) ! value of vapour flux of HDO
160	! REAL h2oflux(ngrid) ! value of vapour flux of H2O
161	REAL old_h2o_vap(ngrid) ! traceur d'eau avant traitement
162	REAL saved_h2o_vap(ngrid) ! traceur d'eau avant traitement
163
164	REAL kmixmin
165
166	c Argument added for surface water ice budget:
167	REAL,INTENT(IN) :: watercap(ngrid,nslope)
168	REAL,INTENT(OUT) :: dwatercap_dif(ngrid,nslope)
169
170	c Subtimestep to compute h2o latent heat flux:
171	REAL :: dtmax = 0.5 ! subtimestep temp criterion
172	INTEGER tsub ! adaptative subtimestep (seconds)
173	REAL subtimestep !ptimestep/nsubtimestep
174	INTEGER nsubtimestep(ngrid) ! number of subtimestep (int)
175
176	c Mass-variation scheme :
177	c ~~~~~~~
178
179	INTEGER j,l
180	REAL zcondicea(ngrid,nlay)
181	REAL zt(ngrid,nlay),ztcond(ngrid,nlay+1)
182	REAL betam(ngrid,nlay),dmice(ngrid,nlay)
183	REAL pdtc(ngrid,nlay)
184	REAL zhs(ngrid,nlay)
185	REAL,SAVE :: ccond
186
187	!$OMP THREADPRIVATE(ccond)
188
189	c Theta_m formulation for mass-variation scheme :
190	c ~~~~~~~
191
192	INTEGER,SAVE :: ico2
193	INTEGER llnt(ngrid)
194	REAL,SAVE :: m_co2, m_noco2, A , B
195	REAL vmr_co2(ngrid,nlay)
196	REAL qco2,mmean
197
198	!$OMP THREADPRIVATE(ico2,m_co2,m_noco2,A,B)
199
200	REAL,INTENT(OUT) :: sensibFlux(ngrid)
201
202	!!MARGAUX
203	REAL DoH_vap(ngrid,nlay)
204	!! Sub-grid scale slopes
205	REAL :: pdqsdif_tmp(ngrid,nq) ! Temporary for dust lifting
206	REAL :: watercap_tmp(ngrid)
207	REAL :: zq_slope_vap(ngrid,nlay,nq,nslope)
208	REAL :: zq_tmp_vap(ngrid,nlay,nq)
209	REAL :: ptsrf_tmp(ngrid)
210	REAL :: pqsurf_tmp(ngrid,nq)
211	REAL :: pdqsdif_tmphdo(ngrid,nq)
212	REAL :: qsat_tmp(ngrid)
213	INTEGER :: indmax
214
215	character*2 str2
216
217	c ** un petit test de coherence
218	c --------------------------
219
220	! AS: OK firstcall absolute
221	IF (firstcall) THEN
222	c To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
223	bcond=1./tcond1mb
224	acond=r/latcond
225	ccond=cpp/(g*latcond)
226	PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
227	PRINT*,' acond,bcond,ccond',acond,bcond,ccond
228
229
230	ico2=0
231
232	c Prepare Special treatment if one of the tracer is CO2 gas
233	do iq=1,nq
234	if (noms(iq).eq."co2") then
235	ico2=iq
236	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
237	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
238	c Compute A and B coefficient use to compute
239	c mean molecular mass Mair defined by
240	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
241	c 1/Mair = A*q(ico2) + B
242	A =(1/m_co2 - 1/m_noco2)
243	B=1/m_noco2
244	endif
245	enddo
246
247	firstcall=.false.
248	ENDIF
249
250	DO ig = 1,ngrid
251	ptsrf_tmp(ig) = ptsrf(ig,major_slope(ig))
252	pqsurf_tmp(ig,:) = pqsurf(ig,:,major_slope(ig))
253	ENDDO
254
255	c-----------------------------------------------------------------------
256	c 1. initialisation
257	c -----------------
258
259	nlev=nlay+1
260
261	! initialize output tendencies to zero:
262	pdudif(1:ngrid,1:nlay)=0
263	pdvdif(1:ngrid,1:nlay)=0
264	pdhdif(1:ngrid,1:nlay)=0
265	pdtsrf(1:ngrid,1:nslope)=0
266	zdtsrf(1:ngrid,1:nslope)=0
267	surf_h2o_lh(1:ngrid,1:nslope)=0
268	zsurf_h2o_lh(1:ngrid,1:nslope)=0
269	pdqdif(1:ngrid,1:nlay,1:nq)=0
270	pdqsdif(1:ngrid,1:nq,1:nslope)=0
271	pdqsdif_tmp(1:ngrid,1:nq)=0
272	zdqsdif(1:ngrid)=0
273	dwatercap_dif(1:ngrid,1:nslope)=0
274
275	c ** calcul de rhodz et dtrho/dz=dtrho*2 g/dp
276	c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
277	c ----------------------------------------
278
279	DO ilay=1,nlay
280	DO ig=1,ngrid
281	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
282	! Mass variation scheme:
283	betam(ig,ilay)=-za(ig,ilay)latcond/(cppppopsk(ig,ilay))
284	ENDDO
285	ENDDO
286
287	zcst1=4.gptimestep/(r*r)
288	DO ilev=2,nlev-1
289	DO ig=1,ngrid
290	zb0(ig,ilev)=pplev(ig,ilev)*
291	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
292	s (ph(ig,ilev-1)+ph(ig,ilev))
293	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
294	s (pplay(ig,ilev-1)-pplay(ig,ilev))
295	ENDDO
296	ENDDO
297	DO ig=1,ngrid
298	zb0(ig,1)=ptimesteppplev(ig,1)/(rptsrf_tmp(ig))
299	ENDDO
300
301	c ** diagnostique pour l'initialisation
302	c ----------------------------------
303
304	IF(lecrit) THEN
305	ig=ngrid/2+1
306	PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
307	DO ilay=1,nlay
308	WRITE(*,'(6f11.5)')
309	s .01pplay(ig,ilay),.001pzlay(ig,ilay),
310	s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
311	ENDDO
312	PRINT*,'Pression (mbar) ,altitude (km),zb'
313	DO ilev=1,nlay
314	WRITE(*,'(3f15.7)')
315	s .01pplev(ig,ilev),.001pzlev(ig,ilev),
316	s zb0(ig,ilev)
317	ENDDO
318	ENDIF
319
320	c -----------------------------------
321	c Potential Condensation temperature:
322	c -----------------------------------
323
324	c Compute CO2 Volume mixing ratio
325	c -------------------------------
326	if (callcond.and.(ico2.ne.0)) then
327	DO ilev=1,nlay
328	DO ig=1,ngrid
329	qco2=MAX(1.E-30
330	& ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep)
331	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
332	mmean=1/(A*qco2 +B)
333	vmr_co2(ig,ilev) = qco2*mmean/m_co2
334	ENDDO
335	ENDDO
336	else
337	DO ilev=1,nlay
338	DO ig=1,ngrid
339	vmr_co2(ig,ilev)=0.95
340	ENDDO
341	ENDDO
342	end if
343
344	c forecast of atmospheric temperature zt and frost temperature ztcond
345	c --------------------------------------------------------------------
346
347	if (callcond) then
348	DO ilev=1,nlay
349	DO ig=1,ngrid
350	ztcond(ig,ilev)=
351	& 1./(bcond-acondlog(.01vmr_co2(ig,ilev)*pplay(ig,ilev)))
352	if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica
353	! zhcond(ig,ilev) =
354	! & (1./(bcond-acondlog(.0095pplay(ig,ilev))))/ppopsk(ig,ilev)
355	zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev)
356	END DO
357	END DO
358	ztcond(:,nlay+1)=ztcond(:,nlay)
359	else
360	zhcond(:,:) = 0
361	ztcond(:,:) = 0
362	end if
363
364
365	c-----------------------------------------------------------------------
366	c 2. ajout des tendances physiques
367	c -----------------------------
368
369	DO ilev=1,nlay
370	DO ig=1,ngrid
371	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
372	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
373	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
374	! zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
375	ENDDO
376	ENDDO
377
378	zq(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)+
379	& pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
380
381	c-----------------------------------------------------------------------
382	c 3. schema de turbulence
383	c --------------------
384
385	c ** source d'energie cinetique turbulente a la surface
386	c (condition aux limites du schema de diffusion turbulente
387	c dans la couche limite
388	c ---------------------
389
390	CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf_tmp
391	& ,ph,zcdv_true,zcdh_true)
392
393	zu2(:)=pu(:,1)pu(:,1)+pv(:,1)pv(:,1)
394
395	IF (callrichsl) THEN
396	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+
397	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
398	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+
399	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
400
401	ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+
402	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
403
404	tstar(:)=0.
405	DO ig=1,ngrid
406	IF (zcdh_true(ig) .ne. 0.) THEN ! When Cd=Ch=0, u=t=0
407	tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
408	& *zcdh(ig)/ustar(ig)
409	ENDIF
410	ENDDO
411
412	ELSE
413	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance
414	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance
415	ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
416	tstar(:)=(ph(:,1)-ptsrf_tmp(:))
417	& *zcdh_true(:)/sqrt(zcdv_true(:))
418	ENDIF
419
420	! Some usefull diagnostics for the new surface layer parametrization :
421
422	! call write_output('vdifc_zcdv_true',
423	! & 'momentum drag','no units',
424	! & zcdv_true(:))
425	! call write_output('vdifc_zcdh_true',
426	! & 'heat drag','no units',
427	! & zcdh_true(:))
428	! call write_output('vdifc_ust',
429	! & 'friction velocity','m/s',ust(:))
430	! call write_output('vdifc_tst',
431	! & 'friction temperature','K',tst(:))
432	! call write_output('vdifc_zcdv',
433	! & 'aerodyn momentum conductance','m/s',
434	! & zcdv(:))
435	! call write_output('vdifc_zcdh',
436	! & 'aerodyn heat conductance','m/s',
437	! & zcdh(:))
438
439	c ** schema de diffusion turbulente dans la couche limite
440	c ----------------------------------------------------
441	IF (.not. callyamada4) THEN
442
443	CALL vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
444	& ,pu,pv,ph,zcdv_true
445	& ,pq2,zkv,zkh,zq)
446
447	ELSE
448
449	pt(:,:)=ph(:,:)*ppopsk(:,:)
450	CALL yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt
451	s ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true,pq2,zkv,zkh,zkq,ustar
452	s ,9)
453	ENDIF
454
455	if ((doubleq).and.(ngrid.eq.1)) then
456	kmixmin = 80. !80.! minimum eddy mix coeff in 1D
457	do ilev=1,nlay
458	do ig=1,ngrid
459	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
460	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
461	end do
462	end do
463	end if
464
465	c ** diagnostique pour le schema de turbulence
466	c -----------------------------------------
467
468	IF(lecrit) THEN
469	PRINT*
470	PRINT*,'Diagnostic for the vertical turbulent mixing'
471	PRINT*,'Cd for momentum and potential temperature'
472
473	PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
474	PRINT*,'Mixing coefficient for momentum and pot.temp.'
475	DO ilev=1,nlay
476	PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
477	ENDDO
478	ENDIF
479
480	c-----------------------------------------------------------------------
481	c 4. inversion pour l'implicite sur u
482	c --------------------------------
483
484	c ** l'equation est
485	c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
486	c avec
487	c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
488	c et
489	c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
490	c donc les entrees sont /zcdv/ pour la condition a la limite sol
491	c et /zkv/ = Ku
492
493	zb(1:ngrid,2:nlay)=zkv(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
494	zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
495
496	DO ig=1,ngrid
497	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
498	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
499	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
500	ENDDO
501
502	DO ilay=nlay-1,1,-1
503	DO ig=1,ngrid
504	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
505	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
506	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
507	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
508	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
509	ENDDO
510	ENDDO
511
512	DO ig=1,ngrid
513	zu(ig,1)=zc(ig,1)
514	ENDDO
515	DO ilay=2,nlay
516	DO ig=1,ngrid
517	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
518	ENDDO
519	ENDDO
520
521	c-----------------------------------------------------------------------
522	c 5. inversion pour l'implicite sur v
523	c --------------------------------
524
525	c ** l'equation est
526	c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
527	c avec
528	c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
529	c et
530	c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
531	c donc les entrees sont /zcdv/ pour la condition a la limite sol
532	c et /zkv/ = Kv
533
534	DO ig=1,ngrid
535	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
536	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
537	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
538	ENDDO
539
540	DO ilay=nlay-1,1,-1
541	DO ig=1,ngrid
542	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
543	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
544	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
545	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
546	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
547	ENDDO
548	ENDDO
549
550	DO ig=1,ngrid
551	zv(ig,1)=zc(ig,1)
552	ENDDO
553	DO ilay=2,nlay
554	DO ig=1,ngrid
555	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
556	ENDDO
557	ENDDO
558
559	c-----------------------------------------------------------------------
560	c 6. inversion pour l'implicite sur h sans oublier le couplage
561	c avec le sol (conduction)
562	c ------------------------
563
564	c ** l'equation est
565	c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
566	c avec
567	c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
568	c et
569	c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
570	c donc les entrees sont /zcdh/ pour la condition de raccord au sol
571	c et /zkh/ = Kh
572	c -------------
573
574	c Mass variation scheme:
575	zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
576	zb(1:ngrid,1)=zcdh(1:ngrid)*zb0(1:ngrid,1)
577
578	c on initialise dm c
579
580	pdtc(:,:)=0.
581	zt(:,:)=0.
582	dmice(:,:)=0.
583
584	c ** calcul de (d Planck / dT) a la temperature d'interface
585	c ------------------------------------------------------
586
587	z4st=4.5.67e-8ptimestep
588	IF (tke_heat_flux .eq. 0.) THEN
589	DO ig=1,ngrid
590	indmax = major_slope(ig)
591	zdplanck(ig)=z4stpemis(ig,indmax)ptsrf(ig,indmax)*
592	& ptsrf(ig,indmax)*ptsrf(ig,indmax)
593	ENDDO
594	ELSE
595	zdplanck(:)=0.
596	ENDIF
597
598	! calcul de zc et zd pour la couche top en prenant en compte le terme
599	! de variation de masse (on fait une boucle pour que \E7a converge)
600
601	! Identification des points de grilles qui ont besoin de la correction
602
603	llnt(:)=1
604	IF (.not.turb_resolved) THEN
605	IF (callcond) THEN
606	DO ig=1,ngrid
607	DO l=1,nlay
608	if(zh(ig,l) .lt. zhcond(ig,l)) then
609	llnt(ig)=300
610	! 200 and 100 do not go beyond month 9 with normal dissipation
611	goto 5
612	endif
613	ENDDO
614	5 continue
615	ENDDO
616	ENDIF
617
618	ENDIF
619
620
621
622	DO ig=1,ngrid
623	indmax = major_slope(ig)
624	! Initialization of z1 and zd, which do not depend on dmice
625
626	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
627	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
628
629	DO ilay=nlay-1,1,-1
630	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
631	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
632	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
633	ENDDO
634
635	! Convergence loop
636
637	DO j=1,llnt(ig)
638
639	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
640	zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
641	& -betam(ig,nlay)*dmice(ig,nlay)
642	zc(ig,nlay)=zc(ig,nlay)*z1(ig)
643	! zd(ig,nlay)=zb(ig,nlay)*z1(ig)
644
645	! calcul de zc et zd pour les couches du haut vers le bas
646
647	DO ilay=nlay-1,1,-1
648	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
649	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
650	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
651	$ zb(ig,ilay+1)*zc(ig,ilay+1)-
652	$ betam(ig,ilay)dmice(ig,ilay))z1(ig)
653	! zd(ig,ilay)=zb(ig,ilay)*z1(ig)
654	ENDDO
655
656	c ** calcul de la temperature_d'interface et de sa tendance.
657	c on ecrit que la somme des flux est nulle a l'interface
658	c a t + \delta t,
659	c c'est a dire le flux radiatif a {t + \delta t}
660	c + le flux turbulent a {t + \delta t}
661	c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
662	c (notation K dt = /cpp*b/)
663	c + le flux dans le sol a t
664	c + l'evolution du flux dans le sol lorsque la temperature d'interface
665	c passe de sa valeur a t a sa valeur a {t + \delta t}.
666	c ----------------------------------------------------
667
668	z1(ig)=pcapcal(ig,indmax)*ptsrf(ig,indmax)
669	s + cppzb(ig,1)zc(ig,1)
670	s + zdplanck(ig)*ptsrf(ig,indmax)
671	s + pfluxsrf(ig,indmax)*ptimestep
672	z2(ig)= pcapcal(ig,indmax)+cppzb(ig,1)(1.-zd(ig,1))
673	s +zdplanck(ig)
674	ztsrf2(ig)=z1(ig)/z2(ig)
675	! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop
676	zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
677
678	c ** et a partir de la temperature au sol on remonte
679	c -----------------------------------------------
680
681	DO ilay=2,nlay
682	zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
683	ENDDO
684	DO ilay=1,nlay
685	zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
686	ENDDO
687
688	c Condensation/sublimation in the atmosphere
689	c ------------------------------------------
690	c (computation of zcondicea and dmice)
691
692	IF (.NOT. co2clouds) then
693	DO l=nlay , 1, -1
694	IF(zt(ig,l).LT.ztcond(ig,l)) THEN
695	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
696	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
697	& ccondpdtc(ig,l)
698	dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
699	END IF
700	ENDDO
701	ELSE
702	DO l=nlay , 1, -1
703	zcondicea(ig,l)= 0.!pcondicea_co2microp(ig,l)*
704	c & (pplev(ig,l) - pplev(ig,l+1))/g
705	dmice(ig,l)= 0.!dmice(ig,l) + zcondicea(ig,l)*ptimestep
706	pdtc(ig,l)=0.
707	ENDDO
708	ENDIF
709
710	ENDDO!of Do j=1,XXX
711	pdtsrf(ig,indmax)=(ztsrf2(ig)-ptsrf(ig,indmax))/ptimestep
712	ENDDO !of Do ig=1,ngrid
713
714	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
715	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zhs(ig,1)-ztsrf2(ig))
716	ENDDO
717
718	c Now implicit sheme on each sub-grid subslope:
719	IF (nslope.ne.1) then
720	DO islope=1,nslope
721	DO ig=1,ngrid
722	IF(islope.ne.major_slope(ig)) then
723	IF (tke_heat_flux .eq. 0.) THEN
724	zdplanck(ig)=z4stpemis(ig,islope)ptsrf(ig,islope)**3
725	ELSE
726	zdplanck(ig) = 0.
727	ENDIF
728	z1(ig)=pcapcal(ig,islope)*ptsrf(ig,islope)
729	s + cppzb(ig,1)zc(ig,1)
730	s + zdplanck(ig)*ptsrf(ig,islope)
731	s + pfluxsrf(ig,islope)*ptimestep
732	z2(ig)= pcapcal(ig,islope)+cppzb(ig,1)(1.-zd(ig,1))
733	s +zdplanck(ig)
734	ztsrf2(ig)=z1(ig)/z2(ig)
735	pdtsrf(ig,islope)=(ztsrf2(ig)-ptsrf(ig,islope))/ptimestep
736	ENDIF ! islope != indmax
737	ENDDO ! ig
738	ENDDO !islope
739	ENDIF !nslope.ne.1
740
741	c-----------------------------------------------------------------------
742	c TRACERS
743	c -------
744
745	c Using the wind modified by friction for lifting and sublimation
746	c ----------------------------------------------------------------
747
748	! This is computed above and takes into account surface-atmosphere flux
749	! enhancement by subgrid gustiness and atmospheric-stability related
750	! variations of transfer coefficients.
751
752	! DO ig=1,ngrid
753	! zu2(ig)=zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
754	! zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
755	! zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
756	! ENDDO
757
758	c Calcul du flux vertical au bas de la premiere couche (dust) :
759	c -----------------------------------------------------------
760	do ig=1,ngrid
761	rho(ig) = zb0(ig,1) /ptimestep
762	c zb(ig,1) = 0.
763	end do
764	c Dust lifting:
765	if (lifting) then
766	#ifndef MESOSCALE
767	if (doubleq.AND.submicron) then
768	do ig=1,ngrid
769	c if(qsurf(ig,igcm_co2).lt.1) then
770	pdqsdif_tmp(ig,igcm_dust_mass) =
771	& -alpha_lift(igcm_dust_mass)
772	pdqsdif_tmp(ig,igcm_dust_number) =
773	& -alpha_lift(igcm_dust_number)
774	pdqsdif_tmp(ig,igcm_dust_submicron) =
775	& -alpha_lift(igcm_dust_submicron)
776	c end if
777	end do
778	else if (doubleq) then
779	if (dustinjection.eq.0) then !injection scheme 0 (old)
780	!or 2 (injection in CL)
781	do ig=1,ngrid
782	if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2
783	pdqsdif_tmp(ig,igcm_dust_mass) =
784	& -alpha_lift(igcm_dust_mass)
785	pdqsdif_tmp(ig,igcm_dust_number) =
786	& -alpha_lift(igcm_dust_number)
787	end if
788	end do
789	elseif(dustinjection.eq.1)then ! dust injection scheme = 1 injection from surface
790	do ig=1,ngrid
791	if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2
792	IF((ti_injection_sol.LE.local_time(ig)).and.
793	& (local_time(ig).LE.tf_injection_sol)) THEN
794	if (rdstorm) then !Rocket dust storm scheme
795	pdqsdif_tmp(ig,igcm_stormdust_mass) =
796	& -alpha_lift(igcm_stormdust_mass)
797	& *dustliftday(ig)
798	pdqsdif_tmp(ig,igcm_stormdust_number) =
799	& -alpha_lift(igcm_stormdust_number)
800	& *dustliftday(ig)
801	pdqsdif_tmp(ig,igcm_dust_mass)= 0.
802	pdqsdif_tmp(ig,igcm_dust_number)= 0.
803	else
804	pdqsdif_tmp(ig,igcm_dust_mass)=
805	& -dustliftday(ig)*
806	& alpha_lift(igcm_dust_mass)
807	pdqsdif_tmp(ig,igcm_dust_number)=
808	& -dustliftday(ig)*
809	& alpha_lift(igcm_dust_number)
810	endif
811	if (submicron) then
812	pdqsdif_tmp(ig,igcm_dust_submicron) = 0.
813	endif ! if (submicron)
814	ELSE ! outside dust injection time frame
815	pdqsdif_tmp(ig,igcm_dust_mass)= 0.
816	pdqsdif_tmp(ig,igcm_dust_number)= 0.
817	if (rdstorm) then
818	pdqsdif_tmp(ig,igcm_stormdust_mass)= 0.
819	pdqsdif_tmp(ig,igcm_stormdust_number)= 0.
820	end if
821	ENDIF
822
823	end if ! of if(qsurf(ig,igcm_co2).lt.1)
824	end do
825	endif ! end if dustinjection
826	else if (submicron) then
827	do ig=1,ngrid
828	pdqsdif_tmp(ig,igcm_dust_submicron) =
829	& -alpha_lift(igcm_dust_submicron)
830	end do
831	else
832	#endif
833	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,
834	& pqsurf_tmp(:,igcm_co2),pdqsdif_tmp)
835	#ifndef MESOSCALE
836	endif !doubleq.AND.submicron
837	#endif
838	else
839	pdqsdif_tmp(1:ngrid,1:nq) = 0.
840	end if
841
842	c OU calcul de la valeur de q a la surface (water) :
843	c ----------------------------------------
844
845	c Inversion pour l'implicite sur q
846	c Cas des traceurs qui ne sont pas h2o_vap
847	c h2o_vap est traite plus loin avec un sous pas de temps
848	c hdo_vap est traite ensuite car dependant de h2o_vap
849	c --------------------------------
850
851	do iq=1,nq !for all tracers except water vapor
852	if ((.not. water).or.(.not. iq.eq.igcm_h2o_vap).or.
853	& (.not. iq.eq.igcm_hdo_vap)) then
854
855
856	zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
857	zb(1:ngrid,1)=0
858
859	DO ig=1,ngrid
860	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
861	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
862	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
863	ENDDO
864
865	DO ilay=nlay-1,2,-1
866	DO ig=1,ngrid
867	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
868	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
869	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
870	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
871	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
872	ENDDO
873	ENDDO
874
875	if ((iq.eq.igcm_h2o_ice)
876	$ .or. (hdo.and.(iq.eq.igcm_hdo_ice) )) then
877
878	DO ig=1,ngrid
879	z1(ig)=1./(za(ig,1)+zb(ig,1)+
880	$ zb(ig,2)*(1.-zd(ig,2)))
881	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
882	$ zb(ig,2)zc(ig,2)) z1(ig) !special case h2o_ice
883	ENDDO
884	else ! every other tracer
885	DO ig=1,ngrid
886	z1(ig)=1./(za(ig,1)+zb(ig,1)+
887	$ zb(ig,2)*(1.-zd(ig,2)))
888	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
889	$ zb(ig,2)*zc(ig,2) +
890	$ (-pdqsdif_tmp(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
891	ENDDO
892	endif !((iq.eq.igcm_h2o_ice)
893	c Starting upward calculations for simple mixing of tracer (dust)
894	DO ig=1,ngrid
895	zq(ig,1,iq)=zc(ig,1)
896	DO ilay=2,nlay
897	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
898	ENDDO
899	ENDDO
900	DO islope = 1,nslope
901	DO ig = 1,ngrid
902	pdqsdif(ig,iq,islope) = pdqsdif_tmp(ig,iq)
903	& * cos(pi*def_slope_mean(islope)/180.)
904	ENDDO
905	ENDDO
906
907	endif! ((.not. water).or.(.not. iq.eq.igcm_h2o_vap)) then
908	enddo ! of do iq=1,nq
909
910	c --------- h2o_vap --------------------------------
911
912
913	c Traitement de la vapeur d'eau h2o_vap
914	c Utilisation d'un sous pas de temps afin
915	c de decrire le flux de chaleur latente
916
917	do iq=1,nq
918	if ((water).and.(iq.eq.igcm_h2o_vap)) then
919
920	DO islope = 1,nslope
921	DO ig=1,ngrid
922	zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/
923	& cos(pi*def_slope_mean(islope)/180.)
924	watercap_tmp(ig) = watercap(ig,islope)/
925	& cos(pi*def_slope_mean(islope)/180.)
926	ENDDO ! ig=1,ngrid
927
928	c make_tsub : sous pas de temps adaptatif
929	c la subroutine est a la fin du fichier
930
931	call make_tsub(ngrid,pdtsrf(:,islope),zqsurf,
932	& ptimestep,dtmax,watercaptag,
933	& nsubtimestep)
934
935	c Calculation for turbulent exchange with the surface (for ice)
936	c initialization of ztsrf, which is surface temperature in
937	c the subtimestep.
938	saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap)
939
940	DO ig=1,ngrid
941	subtimestep = ptimestep/nsubtimestep(ig)
942	ztsrf(ig)=ptsrf(ig,islope) ! +pdtsrf(ig)*subtimestep
943	zq_tmp_vap(ig,:,:) =zq(ig,:,:)
944	c Debut du sous pas de temps
945
946	DO tsub=1,nsubtimestep(ig)
947
948	c C'est parti !
949	zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
950	& /float(nsubtimestep(ig))
951	zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
952	& /float(nsubtimestep(ig))
953	zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1)
954
955	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
956	zc(ig,nlay)=za(ig,nlay)zq_tmp_vap(ig,nlay,iq)z1(ig)
957	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
958	DO ilay=nlay-1,2,-1
959	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
960	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
961	zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+
962	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
963	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
964	ENDDO
965	z1(ig)=1./(za(ig,1)+zb(ig,1)+
966	$ zb(ig,2)*(1.-zd(ig,2)))
967	zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
968	$ zb(ig,2)zc(ig,2)) z1(ig)
969
970	call watersat(1,ztsrf(ig),pplev(ig,1),qsat(ig))
971	old_h2o_vap(ig)=zq_tmp_vap(ig,1,igcm_h2o_vap)
972	zd(ig,1)=zb(ig,1)*z1(ig)
973	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
974
975	zdqsdif(ig)=rho(ig)dryness(ig)zcdv(ig)
976	& *(zq1temp(ig)-qsat(ig))
977	c write(,)'subliming more than available frost: qsurf!'
978
979	if(.not.watercaptag(ig)) then
980	if ((-zdqsdif(ig)*subtimestep)
981	& .gt.(zqsurf(ig))) then
982	c pdqsdif > 0 : ice condensing
983	c pdqsdif < 0 : ice subliming
984	c write(,) "subliming more than available frost: qsurf!"
985	zdqsdif(ig)=
986	& -zqsurf(ig)/subtimestep
987	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
988	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
989	zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,igcm_h2o_vap)+
990	$ zb(ig,2)*zc(ig,2) +
991	$ (-zdqsdif(ig)) subtimestep) z1(ig)
992	zq1temp(ig)=zc(ig,1)
993	endif !if .not.watercaptag(ig)
994	endif ! if sublim more than surface
995
996	c Starting upward calculations for water :
997	c Actualisation de h2o_vap dans le premier niveau
998	zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
999
1000	c Take into account the H2O latent heat impact on the surface temperature
1001	if (latentheat_surfwater) then
1002	lh=(2834.3-0.28*(ztsrf(ig)-To)-
1003	& 0.004(ztsrf(ig)-To)(ztsrf(ig)-To))*1.e+3
1004	zdtsrf(ig,islope)= zdqsdif(ig)*lh /pcapcal(ig,islope)
1005	endif ! (latentheat_surfwater) then
1006
1007	DO ilay=2,nlay
1008	zq_tmp_vap(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)
1009	& *zq_tmp_vap(ig,ilay-1,iq)
1010	ENDDO
1011
1012	c Subtimestep water budget :
1013
1014	ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig,islope)
1015	& + zdtsrf(ig,islope))*subtimestep
1016	zqsurf(ig)= zqsurf(ig)+(
1017	& zdqsdif(ig))*subtimestep
1018
1019	c Monitoring instantaneous latent heat flux in W.m-2 :
1020	zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
1021	& (zdtsrf(ig,islope)*pcapcal(ig,islope))
1022	& *subtimestep
1023
1024	c We ensure that surface temperature can't rise above the solid domain if there
1025	c is still ice on the surface (oldschool)
1026	if(zqsurf(ig)
1027	& +zdqsdif(ig)*subtimestep
1028	& .gt.frost_albedo_threshold) then ! if there is still ice, T cannot exceed To
1029	zdtsrf(ig,islope) = min(zdtsrf(ig,islope),
1030	& (To-ztsrf(ig))/subtimestep) ! ice melt case
1031	endif
1032
1033	c Fin du sous pas de temps
1034	ENDDO ! tsub=1,nsubtimestep
1035
1036	c Integration of subtimestep temp and water budget :
1037	c (btw could also compute the post timestep temp and ice
1038	c by simply adding the subtimestep trend instead of this)
1039	surf_h2o_lh(ig,islope)= zsurf_h2o_lh(ig,islope)/ptimestep
1040	pdtsrf(ig,islope)= (ztsrf(ig) -
1041	& ptsrf(ig,islope))/ptimestep
1042	pdqsdif(ig,igcm_h2o_ice,islope)=
1043	& (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice,islope)/
1044	& cos(pi*def_slope_mean(islope)/180.))
1045	& /ptimestep
1046
1047	c if subliming more than qsurf(ice) and on watercaptag, water
1048	c sublimates from watercap reservoir (dwatercap_dif is <0)
1049	if(watercaptag(ig)) then
1050	if ((-pdqsdif(ig,igcm_h2o_ice,islope)*ptimestep)
1051	& .gt.(pqsurf(ig,igcm_h2o_ice,islope)
1052	& /cos(pi*def_slope_mean(islope)/180.))) then
1053	dwatercap_dif(ig,islope)=
1054	& pdqsdif(ig,igcm_h2o_ice,islope)+
1055	& (pqsurf(ig,igcm_h2o_ice,islope) /
1056	& cos(pi*def_slope_mean(islope)/180.))/ptimestep
1057	pdqsdif(ig,igcm_h2o_ice,islope)=
1058	& - (pqsurf(ig,igcm_h2o_ice,islope)/
1059	& cos(pi*def_slope_mean(islope)/180.))/ptimestep
1060	endif! ((-pdqsdif(ig)*ptimestep)
1061	endif !(watercaptag(ig)) then
1062	zq_slope_vap(ig,:,:,islope) = zq_tmp_vap(ig,:,:)
1063	ENDDO ! of DO ig=1,ngrid
1064	ENDDO ! islope
1065
1066	c Some grid box averages: interface with the atmosphere
1067	DO ig = 1,ngrid
1068	DO ilay = 1,nlay
1069	zq(ig,ilay,iq) = 0.
1070	DO islope = 1,nslope
1071	zq(ig,ilay,iq) = zq(ig,ilay,iq) +
1072	$ zq_slope_vap(ig,ilay,iq,islope) *
1073	$ subslope_dist(ig,islope)
1074	ENDDO
1075	ENDDO
1076	ENDDO
1077
1078	! Recompute values in kg/m^2 slopped
1079	DO ig = 1,ngrid
1080	DO islope = 1,nslope
1081	pdqsdif(ig,igcm_h2o_ice,islope) =
1082	& pdqsdif(ig,igcm_h2o_ice,islope)
1083	& * cos(pi*def_slope_mean(islope)/180.)
1084
1085	dwatercap_dif(ig,islope) =
1086	& dwatercap_dif(ig,islope)
1087	& * cos(pi*def_slope_mean(islope)/180.)
1088	ENDDO
1089	ENDDO
1090
1091	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
1092
1093	c --------- end of h2o_vap ----------------------------
1094
1095	c --------- hdo_vap -----------------------------------
1096
1097	c hdo_ice has already been with along h2o_ice
1098	c amongst "normal" tracers (ie not h2o_vap)
1099
1100	if (hdo.and.(iq.eq.igcm_hdo_vap)) then
1101	zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
1102	zb(1:ngrid,1)=0
1103
1104	DO ig=1,ngrid
1105	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
1106	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
1107	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
1108	ENDDO
1109
1110	DO ilay=nlay-1,2,-1
1111	DO ig=1,ngrid
1112	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
1113	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
1114	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
1115	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
1116	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
1117	ENDDO
1118	ENDDO
1119	hdoflux_meshavg(:) = 0.
1120	DO islope = 1,nslope
1121
1122	pdqsdif_tmphdo(:,:) = pdqsdif(:,:,islope)
1123	& /cos(pi*def_slope_mean(islope)/180.)
1124
1125	call watersat(ngrid,pdtsrf(:,islope)*ptimestep +
1126	& ptsrf(:,islope),pplev(:,1),qsat_tmp)
1127
1128	CALL hdo_surfex(ngrid,nlay,nq,ptimestep,
1129	& zt,pplay,zq,pqsurf(:,:,islope),
1130	& saved_h2o_vap,qsat_tmp,
1131	& pdqsdif_tmphdo,
1132	& dwatercap_dif(:,islope)/cos(pi*def_slope_mean(islope)/180.),
1133	& hdoflux(:,islope))
1134
1135	pdqsdif(:,:,islope) = pdqsdif_tmphdo(:,:) *
1136	& cos(pi*def_slope_mean(islope)/180.)
1137	DO ig = 1,ngrid
1138	hdoflux_meshavg(ig) = hdoflux_meshavg(ig) +
1139	& hdoflux(ig,islope)*subslope_dist(ig,islope)
1140
1141	ENDDO !ig = 1,ngrid
1142	ENDDO !islope = 1,nslope
1143
1144	DO ig=1,ngrid
1145	z1(ig)=1./(za(ig,1)+zb(ig,1)+
1146	$ zb(ig,2)*(1.-zd(ig,2)))
1147	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
1148	$ zb(ig,2)*zc(ig,2) +
1149	$ (-hdoflux_meshavg(ig)) ptimestep) z1(ig) !tracer flux from surface
1150	ENDDO
1151
1152	DO ig=1,ngrid
1153	zq(ig,1,iq)=zc(ig,1)
1154	DO ilay=2,nlay
1155	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
1156	ENDDO
1157	ENDDO
1158	endif ! (hdo.and.(iq.eq.igcm_hdo_vap))
1159
1160	c --------- end of hdo ----------------------------
1161
1162	enddo ! of do iq=1,nq
1163
1164	c --------- end of tracers ----------------------------
1165
1166	call write_output("surf_h2o_lh",
1167	& "Ground ice latent heat flux",
1168	& "W.m-2",surf_h2o_lh(:,iflat))
1169
1170	C Diagnostic output for HDO
1171	! if (hdo) then
1172	! CALL write_output('hdoflux',
1173	! & 'hdoflux',
1174	! & ' ',hdoflux_meshavg(:))
1175	! CALL write_output('h2oflux',
1176	! & 'h2oflux',
1177	! & ' ',h2oflux(:))
1178	! endif
1179
1180	c-----------------------------------------------------------------------
1181	c 8. calcul final des tendances de la diffusion verticale
1182	c ----------------------------------------------------
1183
1184	DO ilev = 1, nlay
1185	DO ig=1,ngrid
1186	pdudif(ig,ilev)=( zu(ig,ilev)-
1187	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
1188	pdvdif(ig,ilev)=( zv(ig,ilev)-
1189	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
1190	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
1191	$ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
1192	pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
1193	ENDDO
1194	ENDDO
1195
1196	pdqdif(1:ngrid,1:nlay,1:nq)=(zq(1:ngrid,1:nlay,1:nq)-
1197	& (pq(1:ngrid,1:nlay,1:nq)
1198	& +pdqfi(1:ngrid,1:nlay,1:nq)
1199	& *ptimestep))/ptimestep
1200
1201	c ** diagnostique final
1202	c ------------------
1203
1204	IF(lecrit) THEN
1205	PRINT*,'In vdif'
1206	PRINT*,'Ts (t) and Ts (t+st)'
1207	WRITE(*,'(a10,3a15)')
1208	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
1209	PRINT*,ptsrf(ngrid/2+1,:),ztsrf2(ngrid/2+1)
1210	DO ilev=1,nlay
1211	WRITE(*,'(4f15.7)')
1212	s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
1213	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
1214
1215	ENDDO
1216	ENDIF
1217
1218	END SUBROUTINE vdifc
1219
1220	c====================================
1221
1222	SUBROUTINE make_tsub(naersize,dtsurf,qsurf,ptimestep,
1223	$ dtmax,watercaptag,ntsub)
1224
1225	c Pas de temps adaptatif en estimant le taux de sublimation
1226	c et en adaptant avec un critere "dtmax" du chauffage a accomoder
1227	c dtmax est regle empiriquement (pour l'instant) a 0.5 K
1228
1229	integer,intent(in) :: naersize
1230	real,intent(in) :: dtsurf(naersize)
1231	real,intent(in) :: qsurf(naersize)
1232	logical,intent(in) :: watercaptag(naersize)
1233	real,intent(in) :: ptimestep
1234	real,intent(in) :: dtmax
1235	real :: ztsub(naersize)
1236	integer :: i
1237	integer,intent(out) :: ntsub(naersize)
1238
1239	do i=1,naersize
1240	if ((qsurf(i).eq.0).and.
1241	& (.not.watercaptag(i))) then
1242	ntsub(i) = 1
1243	else
1244	ztsub(i) = ptimestep * dtsurf(i) / dtmax
1245	ntsub(i) = ceiling(abs(ztsub(i)))
1246	endif ! (qsurf(i).eq.0) then
1247	c
1248	c write(78,), dtsurfptimestep, dtsurf, ntsub
1249	enddo! 1=1,ngrid
1250
1251
1252
1253	END SUBROUTINE make_tsub
1254	END MODULE vdifc_mod

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

Download in other formats:

Original Format