Context Navigation

vdifc_mod.F @ 2156

Last change on this file since 2156 was 2080, checked in by mvals, 6 years ago

Mars GCM:

typo in vdifc.F (visible in debug mode) in the case of rdstorm=false, condition for defining pdqsdif(stormdust_mass/number)
typo in physiq_mod.F (visible in debug mode): wrong indices for tracers in aerosol(:,:,:)

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

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

Download in other formats:

Original Format