Context Navigation

vdifc_mod.F @ 2541

Last change on this file since 2541 was 2531, checked in by emillour, 4 years ago
Mars GCM: Bug (wrong usage of surface pressure) fix in vdifc (bug was introduced in r2515) EM
File size: 38.4 KB

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

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

Download in other formats:

Original Format