Context Navigation

vdifc_mod.F @ 2800

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

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

Download in other formats:

Original Format