Context Navigation

vdifc_mod.F @ 2393

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

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

Download in other formats:

Original Format