Context Navigation

vdifc_mod.F @ 2079

Last change on this file since 2079 was 1996, checked in by emillour, 7 years ago

Mars physics:

Turn watersat into a module.

CO2 cloud updates:

compute co2 condensation tendencies in the co2 cloud scheme and pass them on to vdifc (for tests; they might not be needed) and adapt newcondens.

DB+EM

File size: 30.7 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 ça 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	pdqsdif(ig,igcm_stormdust_mass)= 0.
745	pdqsdif(ig,igcm_stormdust_number)= 0.
746	ENDIF
747
748	end if ! of if(co2ice(ig).lt.1)
749	end do
750	endif ! end if dustinjection
751	else if (submicron) then
752	do ig=1,ngrid
753	pdqsdif(ig,igcm_dust_submicron) =
754	& -alpha_lift(igcm_dust_submicron)
755	end do
756	else
757	#endif
758	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
759	& pdqsdif)
760	#ifndef MESOSCALE
761	endif !doubleq.AND.submicron
762	#endif
763	else
764	pdqsdif(1:ngrid,1:nq) = 0.
765	end if
766
767	c OU calcul de la valeur de q a la surface (water) :
768	c ----------------------------------------
769	if (water) then
770	call watersat(ngrid,ptsrf,pplev(1,1),qsat)
771	end if
772
773	c Inversion pour l'implicite sur q
774	c --------------------------------
775	do iq=1,nq !for all tracers including stormdust
776	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
777
778	if ((water).and.(iq.eq.igcm_h2o_vap)) then
779	c This line is required to account for turbulent transport
780	c from surface (e.g. ice) to mid-layer of atmosphere:
781	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
782	CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
783	else ! (re)-initialize zb(:,1)
784	zb(1:ngrid,1)=0
785	end if
786
787	DO ig=1,ngrid
788	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
789	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
790	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
791	ENDDO
792
793	DO ilay=nlay-1,2,-1
794	DO ig=1,ngrid
795	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
796	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
797	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
798	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
799	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
800	ENDDO
801	ENDDO
802
803	if (water.and.(iq.eq.igcm_h2o_ice)) then
804	! special case for water ice tracer: do not include
805	! h2o ice tracer from surface (which is set when handling
806	! h2o vapour case (see further down).
807	DO ig=1,ngrid
808	z1(ig)=1./(za(ig,1)+zb(ig,1)+
809	$ zb(ig,2)*(1.-zd(ig,2)))
810	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
811	$ zb(ig,2)zc(ig,2)) z1(ig)
812	ENDDO
813	else ! general case
814	DO ig=1,ngrid
815	z1(ig)=1./(za(ig,1)+zb(ig,1)+
816	$ zb(ig,2)*(1.-zd(ig,2)))
817	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
818	$ zb(ig,2)*zc(ig,2) +
819	$ (-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
820	ENDDO
821	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
822
823	IF ((water).and.(iq.eq.igcm_h2o_vap)) then
824	c Calculation for turbulent exchange with the surface (for ice)
825	DO ig=1,ngrid
826	zd(ig,1)=zb(ig,1)*z1(ig)
827	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
828
829	pdqsdif(ig,igcm_h2o_ice)=rho(ig)dryness(ig)zcdv(ig)
830	& *(zq1temp(ig)-qsat(ig))
831	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
832	END DO
833
834	DO ig=1,ngrid
835	if(.not.watercaptag(ig)) then
836	if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
837	& .gt.pqsurf(ig,igcm_h2o_ice)) then
838	c write(,)'on sublime plus que qsurf!'
839	pdqsdif(ig,igcm_h2o_ice)=
840	& -pqsurf(ig,igcm_h2o_ice)/ptimestep
841	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
842	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
843	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
844	$ zb(ig,2)*zc(ig,2) +
845	$ (-pdqsdif(ig,igcm_h2o_ice)) ptimestep) z1(ig)
846	zq1temp(ig)=zc(ig,1)
847	endif
848	endif ! if (.not.watercaptag(ig))
849	c Starting upward calculations for water :
850	zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
851
852	!c Take into account H2O latent heat in surface energy budget
853	!c We solve dT/dt = (2834.3-0.28(T-To)-0.004(T-To)^2)1e3iceflux/cpp
854	! tsrf_lw(ig) = ptsrf(ig) + pdtsrf(ig) *ptimestep
855	!
856	! alpha = exp(-4abs(T1-T2)pdqsdif(ig,igcm_h2o_ice)
857	! & *ptimestep/pcapcal(ig))
858	!
859	! tsrf_lw(ig) = (tsrf_lw(ig)(T2-alphaT1)+T1T2(alpha-1))
860	! & /(tsrf_lw(ig)(1-alpha)+alphaT2-T1) ! surface temperature at t+1
861	!
862	! pdtsrf(ig) = (tsrf_lw(ig)-ptsrf(ig))/ptimestep
863
864	if(pqsurf(ig,igcm_h2o_ice)
865	& +pdqsdif(ig,igcm_h2o_ice)*ptimestep
866	& .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To
867	& pdtsrf(ig) = min(pdtsrf(ig),(To-ptsrf(ig))/ptimestep) ! ice melt case
868
869	ENDDO ! of DO ig=1,ngrid
870	ELSE
871	c Starting upward calculations for simple mixing of tracer (dust)
872	DO ig=1,ngrid
873	zq(ig,1,iq)=zc(ig,1)
874	ENDDO
875	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
876
877	DO ilay=2,nlay
878	DO ig=1,ngrid
879	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
880	ENDDO
881	ENDDO
882	enddo ! of do iq=1,nq
883	end if ! of if(tracer)
884
885
886	c-----------------------------------------------------------------------
887	c 8. calcul final des tendances de la diffusion verticale
888	c ----------------------------------------------------
889
890	DO ilev = 1, nlay
891	DO ig=1,ngrid
892	pdudif(ig,ilev)=( zu(ig,ilev)-
893	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
894	pdvdif(ig,ilev)=( zv(ig,ilev)-
895	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
896	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
897	$ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
898	pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
899	ENDDO
900	ENDDO
901
902	if (tracer) then
903	DO iq = 1, nq
904	DO ilev = 1, nlay
905	DO ig=1,ngrid
906	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
907	$ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
908	ENDDO
909	ENDDO
910	ENDDO
911	end if
912
913	c ** diagnostique final
914	c ------------------
915
916	IF(lecrit) THEN
917	PRINT*,'In vdif'
918	PRINT*,'Ts (t) and Ts (t+st)'
919	WRITE(*,'(a10,3a15)')
920	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
921	PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
922	DO ilev=1,nlay
923	WRITE(*,'(4f15.7)')
924	s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
925	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
926
927	ENDDO
928	ENDIF
929
930	END SUBROUTINE vdifc
931
932	END MODULE vdifc_mod

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

Download in other formats:

Original Format