Context Navigation

vdifc_mod.F @ 2511

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

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

Download in other formats:

Original Format