Context Navigation

vdifc_mod.F @ 2840

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

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

Download in other formats:

Original Format