Context Navigation

vdifc.F @ 1962

Last change on this file since 1962 was 1944, checked in by emillour, 7 years ago

Mars GCM:
A first step towards 1+1=2 (for now only works without tracers):

store and load "albedo" (surface albedo) and wstar (thermals' max vertical velocity) in physics (re)start file.
turn phyetat0 into a module in the process.

File size: 28.2 KB

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

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

Download in other formats:

Original Format