Context Navigation

vdifc_mod.F @ 2826

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

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

Download in other formats:

Original Format