Context Navigation

vdifc_mod.F @ 2937

Last change on this file since 2937 was 2934, checked in by emillour, 3 years ago
Mars PCM: A first step in cleaning up outputs and adding field definitions for XIOS. EM
File size: 38.9 KB

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

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

Download in other formats:

Original Format