Context Navigation

vdifc_mod.F @ 2270

Last change on this file since 2270 was 2260, checked in by jnaar, 5 years ago

[MARS GCM] Adding a new 2D variable "watercap", perennial ground water ice.
qsurf(ig,igcm_h2o_ice) can no longer be negative.
When watercaptag=true, ie when there is an infinite amount of water avalaible for sublimation,
and there is no more h2o frost, sublimation digs into watercap reservoir.
For now watercap can't grow, even on watercaptag, and has same albedo than qsurf(h2o_ice).
When using old start file, watercap is automatically initiated as 0.
Added watercap in newstart.F aswell.
JN

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

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

Download in other formats:

Original Format