Context Navigation

vdifc_mod.F @ 2218

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

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

Download in other formats:

Original Format