Context Navigation

vdifc_mod.F @ 2613

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

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

Download in other formats:

Original Format