Context Navigation

clmain.F @ 661

Last change on this file since 661 was 661, checked in by Laurent Fairhead, 19 years ago
Bug decouvert par Gerhard, agesno n'est pas mis a jour apres passage dans interface_surf. S'applique aussi aux simulations IPCC ... LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 75.3 KB

Line
1	!
2	! $Header$
3	!
4	c
5	c
6	SUBROUTINE clmain(dtime,itap,date0,pctsrf,pctsrf_new,
7	. t,q,u,v,
8	. jour, rmu0, co2_ppm,
9	. ok_veget, ocean, npas, nexca, ts,
10	. soil_model,cdmmax, cdhmax,
11	. ksta, ksta_ter, ok_kzmin, ftsoil,qsol,
12	cIM BAD . paprs,pplay,radsol,snow,qsurf,evap,albe,alblw,
13	. paprs,pplay,snow,qsurf,evap,albe,alblw,
14	. fluxlat,
15	. rain_f, snow_f, solsw, sollw, sollwdown, fder,
16	. rlon, rlat, cufi, cvfi, rugos,
17	. debut, lafin, agesno,rugoro,
18	. d_t,d_q,d_u,d_v,d_ts,
19	. flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
20	. q2,
21	. dflux_t,dflux_q,
22	. zcoefh,zu1,zv1, t2m, q2m, u10m, v10m,
23	cIM cf. AM : pbl
24	. pblh,capCL,oliqCL,cteiCL,pblT,
25	. therm,trmb1,trmb2,trmb3,plcl,
26	. fqcalving,ffonte, run_off_lic_0,
27	cIM "slab" ocean
28	. flux_o, flux_g, tslab, seaice)
29	cAA . itr, tr, flux_surf, d_tr)
30	cAA REM:
31	cAA-----
32	cAA Tout ce qui a trait au traceurs est dans phytrac maintenant
33	cAA pour l'instant le calcul de la couche limite pour les traceurs
34	cAA se fait avec cltrac et ne tient pas compte de la differentiation
35	cAA des sous-fraction de sol.
36	cAA REM bis :
37	cAA----------
38	cAA Pour pouvoir extraire les coefficient d'echanges et le vent
39	cAA dans la premiere couche, 3 champs supplementaires ont ete crees
40	cAA zcoefh,zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
41	cAA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
42	cAA si les informations des subsurfaces doivent etre prises en compte
43	cAA il faudra sortir ces memes champs en leur ajoutant une dimension,
44	cAA c'est a dire nbsrf (nbre de subsurface).
45	USE ioipsl
46	USE interface_surf
47	IMPLICIT none
48	c======================================================================
49	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
50	c Objet: interface de "couche limite" (diffusion verticale)
51	c Arguments:
52	c dtime----input-R- interval du temps (secondes)
53	c itap-----input-I- numero du pas de temps
54	c date0----input-R- jour initial
55	c t--------input-R- temperature (K)
56	c q--------input-R- vapeur d'eau (kg/kg)
57	c u--------input-R- vitesse u
58	c v--------input-R- vitesse v
59	c ts-------input-R- temperature du sol (en Kelvin)
60	c paprs----input-R- pression a intercouche (Pa)
61	c pplay----input-R- pression au milieu de couche (Pa)
62	c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
63	c rlat-----input-R- latitude en degree
64	c rugos----input-R- longeur de rugosite (en m)
65	c cufi-----input-R- resolution des mailles en x (m)
66	c cvfi-----input-R- resolution des mailles en y (m)
67	c
68	c d_t------output-R- le changement pour "t"
69	c d_q------output-R- le changement pour "q"
70	c d_u------output-R- le changement pour "u"
71	c d_v------output-R- le changement pour "v"
72	c d_ts-----output-R- le changement pour "ts"
73	c flux_t---output-R- flux de chaleur sensible (CpT) J/m2/s (W/m2)
74	c (orientation positive vers le bas)
75	c flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
76	c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
77	c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
78	c dflux_t derive du flux sensible
79	c dflux_q derive du flux latent
80	cIM "slab" ocean
81	c flux_g---output-R- flux glace (pour OCEAN='slab ')
82	c flux_o---output-R- flux ocean (pour OCEAN='slab ')
83	c tslab-in/output-R temperature du slab ocean (en Kelvin) ! uniqmnt pour slab
84	c seaice---output-R- glace de mer (kg/m2) (pour OCEAN='slab ')
85	ccc
86	c ffonte----Flux thermique utilise pour fondre la neige
87	c fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
88	c hauteur de neige, en kg/m2/s
89	cAA on rajoute en output yu1 et yv1 qui sont les vents dans
90	cAA la premiere couche
91	cAA ces 4 variables sont maintenant traites dans phytrac
92	c itr--------input-I- nombre de traceurs
93	c tr---------input-R- q. de traceurs
94	c flux_surf--input-R- flux de traceurs a la surface
95	c d_tr-------output-R tendance de traceurs
96	cIM cf. AM : PBL
97	c trmb1-------deep_cape
98	c trmb2--------inhibition
99	c trmb3-------Point Omega
100	c Cape(klon)-------Cape du thermique
101	c EauLiq(klon)-------Eau liqu integr du thermique
102	c ctei(klon)-------Critere d'instab d'entrainmt des nuages de CL
103	c lcl------- Niveau de condensation
104	c pblh------- HCL
105	c pblT------- T au nveau HCL
106	c======================================================================
107	#include "dimensions.h"
108	#include "dimphy.h"
109	#include "indicesol.h"
110	c$$$ PB ajout pour soil
111	#include "dimsoil.h"
112	#include "iniprint.h"
113	#include "compbl.h"
114	c
115	REAL dtime
116	real date0
117	integer itap
118	REAL t(klon,klev), q(klon,klev)
119	REAL u(klon,klev), v(klon,klev)
120	cIM 230604 BAD REAL radsol(klon) ???
121	REAL paprs(klon,klev+1), pplay(klon,klev)
122	REAL rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
123	REAL d_t(klon, klev), d_q(klon, klev)
124	REAL d_u(klon, klev), d_v(klon, klev)
125	REAL flux_t(klon,klev, nbsrf), flux_q(klon,klev, nbsrf)
126	REAL dflux_t(klon), dflux_q(klon)
127	cIM "slab" ocean
128	REAL flux_o(klon), flux_g(klon)
129	REAL y_flux_o(klon), y_flux_g(klon)
130	REAL tslab(klon), ytslab(klon)
131	REAL seaice(klon), y_seaice(klon)
132	cIM cf JLD
133	REAL y_fqcalving(klon), y_ffonte(klon)
134	REAL fqcalving(klon,nbsrf), ffonte(klon,nbsrf)
135	REAL run_off_lic_0(klon), y_run_off_lic_0(klon)
136
137	REAL flux_u(klon,klev, nbsrf), flux_v(klon,klev, nbsrf)
138	REAL rugmer(klon), agesno(klon,nbsrf),rugoro(klon)
139	REAL cdragh(klon), cdragm(klon)
140	integer jour ! jour de l'annee en cours
141	real rmu0(klon) ! cosinus de l'angle solaire zenithal
142	REAL co2_ppm ! taux CO2 atmosphere
143	LOGICAL debut, lafin, ok_veget
144	character*6 ocean
145	integer npas, nexca
146	cAA INTEGER itr
147	cAA REAL tr(klon,klev,nbtr)
148	cAA REAL d_tr(klon,klev,nbtr)
149	cAA REAL flux_surf(klon,nbtr)
150	c
151	REAL pctsrf(klon,nbsrf)
152	REAL ts(klon,nbsrf)
153	REAL d_ts(klon,nbsrf)
154	REAL snow(klon,nbsrf)
155	REAL qsurf(klon,nbsrf)
156	REAL evap(klon,nbsrf)
157	REAL albe(klon,nbsrf)
158	REAL alblw(klon,nbsrf)
159	c$$$ PB
160	REAL fluxlat(klon,nbsrf)
161	C
162	real rain_f(klon), snow_f(klon)
163	REAL fder(klon)
164	cIM cf. JLD REAL sollw(klon), solsw(klon), sollwdown(klon)
165	REAL sollw(klon,nbsrf), solsw(klon,nbsrf), sollwdown(klon)
166	REAL rugos(klon,nbsrf)
167	C la nouvelle repartition des surfaces sortie de l'interface
168	REAL pctsrf_new(klon,nbsrf)
169	cAA
170	REAL zcoefh(klon,klev)
171	REAL zu1(klon)
172	REAL zv1(klon)
173	cAA
174	c$$$ PB ajout pour soil
175	LOGICAL soil_model
176	cIM ajout seuils cdrm, cdrh
177	REAL cdmmax, cdhmax
178	cIM: 261103
179	REAL ksta, ksta_ter
180	LOGICAL ok_kzmin
181	cIM: 261103
182	REAL ftsoil(klon,nsoilmx,nbsrf)
183	REAL ytsoil(klon,nsoilmx)
184	REAL qsol(klon)
185	c======================================================================
186	EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
187	c======================================================================
188	REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
189	REAL yalb(klon)
190	REAL yalblw(klon)
191	REAL yu1(klon), yv1(klon)
192	real ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
193	real yrain_f(klon), ysnow_f(klon)
194	real ysollw(klon), ysolsw(klon), ysollwdown(klon)
195	real yfder(klon), ytaux(klon), ytauy(klon)
196	REAL yrugm(klon), yrads(klon),yrugoro(klon)
197	c$$$ PB
198	REAL yfluxlat(klon)
199	C
200	REAL y_d_ts(klon)
201	REAL y_d_t(klon, klev), y_d_q(klon, klev)
202	REAL y_d_u(klon, klev), y_d_v(klon, klev)
203	REAL y_flux_t(klon,klev), y_flux_q(klon,klev)
204	REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
205	REAL y_dflux_t(klon), y_dflux_q(klon)
206	REAL ycoefh(klon,klev), ycoefm(klon,klev)
207	REAL yu(klon,klev), yv(klon,klev)
208	REAL yt(klon,klev), yq(klon,klev)
209	REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
210	cAA REAL ytr(klon,klev,nbtr)
211	cAA REAL y_d_tr(klon,klev,nbtr)
212	cAA REAL yflxsrf(klon,nbtr)
213	c
214	c
215	LOGICAL ok_nonloc
216	PARAMETER (ok_nonloc=.FALSE.)
217	REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
218
219	cIM 081204 hcl_Anne ? BEG
220	real yzlay(klon,klev),yzlev(klon,klev+1),yteta(klon,klev)
221	real ykmm(klon,klev+1),ykmn(klon,klev+1)
222	real ykmq(klon,klev+1)
223	real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)
224	real q2diag(klon,klev+1)
225	cIM 081204 real yustar(klon),y_cd_m(klon),y_cd_h(klon)
226	cIM 081204 hcl_Anne ? END
227	c
228	#include "YOMCST.h"
229	#include "YOETHF.h"
230	#include "FCTTRE.h"
231	REAL u1lay(klon), v1lay(klon)
232	REAL delp(klon,klev)
233	INTEGER i, k, nsrf
234	cAA INTEGER it
235	INTEGER ni(klon), knon, j
236	c Introduction d'une variable "pourcentage potentiel" pour tenir compte
237	c des eventuelles apparitions et/ou disparitions de la glace de mer
238	REAL pctsrf_pot(klon,nbsrf)
239
240	c======================================================================
241	REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
242	c======================================================================
243	c
244	c maf pour sorties IOISPL en cas de debugagage
245	c
246	CHARACTER*80 cldebug
247	SAVE cldebug
248	CHARACTER*8 cl_surf(nbsrf)
249	SAVE cl_surf
250	INTEGER nhoridbg, nidbg
251	SAVE nhoridbg, nidbg
252	INTEGER ndexbg(iim*(jjm+1))
253	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian
254	REAL tabindx(klon)
255	REAL debugtab(iim,jjm+1)
256	LOGICAL first_appel
257	SAVE first_appel
258	DATA first_appel/.true./
259	LOGICAL debugindex
260	SAVE debugindex
261	DATA debugindex/.false./
262	integer idayref
263	#include "temps.h"
264	REAL t2m(klon,nbsrf), q2m(klon,nbsrf)
265	REAL u10m(klon,nbsrf), v10m(klon,nbsrf)
266	c
267	REAL yt2m(klon), yq2m(klon), yu10m(klon)
268	REAL yustar(klon)
269	c -- LOOP
270	REAL yu10mx(klon)
271	REAL yu10my(klon)
272	REAL ywindsp(klon)
273	c -- LOOP
274	c
275	REAL yt10m(klon), yq10m(klon)
276	cIM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds physic ce qui
277	c permet de sortir les grdeurs par sous surface)
278	REAL pblh(klon,nbsrf)
279	REAL plcl(klon,nbsrf)
280	REAL capCL(klon,nbsrf)
281	REAL oliqCL(klon,nbsrf)
282	REAL cteiCL(klon,nbsrf)
283	REAL pblT(klon,nbsrf)
284	REAL therm(klon,nbsrf)
285	REAL trmb1(klon,nbsrf)
286	REAL trmb2(klon,nbsrf)
287	REAL trmb3(klon,nbsrf)
288	REAL ypblh(klon)
289	REAL ylcl(klon)
290	REAL ycapCL(klon)
291	REAL yoliqCL(klon)
292	REAL ycteiCL(klon)
293	REAL ypblT(klon)
294	REAL ytherm(klon)
295	REAL ytrmb1(klon)
296	REAL ytrmb2(klon)
297	REAL ytrmb3(klon)
298	REAL y_cd_h(klon), y_cd_m(klon)
299	c REAL ygamt(klon,2:klev) ! contre-gradient pour temperature
300	c REAL ygamq(klon,2:klev) ! contre-gradient pour humidite
301	REAL uzon(klon), vmer(klon)
302	REAL tair1(klon), qair1(klon), tairsol(klon)
303	REAL psfce(klon), patm(klon)
304	c
305	REAL qairsol(klon), zgeo1(klon)
306	REAL rugo1(klon)
307	c
308	LOGICAL zxli ! utiliser un jeu de fonctions simples
309	PARAMETER (zxli=.FALSE.)
310	c
311	REAL zt, zqs, zdelta, zcor
312	REAL t_coup
313	PARAMETER(t_coup=273.15)
314	C
315	character (len = 20) :: modname = 'clmain'
316	LOGICAL check
317	PARAMETER (check=.false.)
318	C
319	if (check) THEN
320	write(,) modname,' klon=',klon
321	call flush(6)
322	endif
323	IF (debugindex .and. first_appel) THEN
324	first_appel=.false.
325	!
326	! initialisation sorties netcdf
327	!
328	idayref = day_ini
329	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
330	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
331	DO i = 1, iim
332	zx_lon(i,1) = rlon(i+1)
333	zx_lon(i,jjm+1) = rlon(i+1)
334	ENDDO
335	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
336	cldebug='sous_index'
337	CALL histbeg(cldebug, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),
338	$ 1,iim,1,jjm
339	$ +1, itau_phy,zjulian,dtime,nhoridbg,nidbg)
340	! no vertical axis
341	cl_surf(1)='ter'
342	cl_surf(2)='lic'
343	cl_surf(3)='oce'
344	cl_surf(4)='sic'
345	DO nsrf=1,nbsrf
346	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim,
347	$ jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
348	END DO
349	CALL histend(nidbg)
350	CALL histsync(nidbg)
351	ENDIF
352
353	DO k = 1, klev ! epaisseur de couche
354	DO i = 1, klon
355	delp(i,k) = paprs(i,k)-paprs(i,k+1)
356	ENDDO
357	ENDDO
358	DO i = 1, klon ! vent de la premiere couche
359	ccc zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
360	zx_alf1 = 1.0
361	zx_alf2 = 1.0 - zx_alf1
362	u1lay(i) = u(i,1)zx_alf1 + u(i,2)zx_alf2
363	v1lay(i) = v(i,1)zx_alf1 + v(i,2)zx_alf2
364	ENDDO
365	c
366	c initialisation:
367	c
368	DO i = 1, klon
369	rugmer(i) = 0.0
370	cdragh(i) = 0.0
371	cdragm(i) = 0.0
372	dflux_t(i) = 0.0
373	dflux_q(i) = 0.0
374	zu1(i) = 0.0
375	zv1(i) = 0.0
376	ENDDO
377	ypct = 0.0
378	yts = 0.0
379	ysnow = 0.0
380	yqsurf = 0.0
381	yalb = 0.0
382	yalblw = 0.0
383	yrain_f = 0.0
384	ysnow_f = 0.0
385	yfder = 0.0
386	ytaux = 0.0
387	ytauy = 0.0
388	ysolsw = 0.0
389	ysollw = 0.0
390	ysollwdown = 0.0
391	yrugos = 0.0
392	yu1 = 0.0
393	yv1 = 0.0
394	yrads = 0.0
395	ypaprs = 0.0
396	ypplay = 0.0
397	ydelp = 0.0
398	yu = 0.0
399	yv = 0.0
400	yt = 0.0
401	yq = 0.0
402	pctsrf_new = 0.0
403	y_flux_u = 0.0
404	y_flux_v = 0.0
405	C$$ PB
406	y_dflux_t = 0.0
407	y_dflux_q = 0.0
408	ytsoil = 999999.
409	yrugoro = 0.
410	c -- LOOP
411	yu10mx = 0.0
412	yu10my = 0.0
413	ywindsp = 0.0
414	c -- LOOP
415	DO nsrf = 1, nbsrf
416	DO i = 1, klon
417	d_ts(i,nsrf) = 0.0
418	ENDDO
419	END DO
420	C§§§ PB
421	yfluxlat=0.
422	flux_t = 0.
423	flux_q = 0.
424	flux_u = 0.
425	flux_v = 0.
426	DO k = 1, klev
427	DO i = 1, klon
428	d_t(i,k) = 0.0
429	d_q(i,k) = 0.0
430	c$$$ flux_t(i,k) = 0.0
431	c$$$ flux_q(i,k) = 0.0
432	d_u(i,k) = 0.0
433	d_v(i,k) = 0.0
434	c$$$ flux_u(i,k) = 0.0
435	c$$$ flux_v(i,k) = 0.0
436	zcoefh(i,k) = 0.0
437	ENDDO
438	ENDDO
439	cAA IF (itr.GE.1) THEN
440	cAA DO it = 1, itr
441	cAA DO k = 1, klev
442	cAA DO i = 1, klon
443	cAA d_tr(i,k,it) = 0.0
444	cAA ENDDO
445	cAA ENDDO
446	cAA ENDDO
447	cAA ENDIF
448
449	c
450	c Boucler sur toutes les sous-fractions du sol:
451	c
452	C Initialisation des "pourcentages potentiels". On considere ici qu'on
453	C peut avoir potentiellementdela glace sur tout le domaine oceanique
454	C (a affiner)
455
456	pctsrf_pot = pctsrf
457	pctsrf_pot(:,is_oce) = 1. - zmasq(:)
458	pctsrf_pot(:,is_sic) = 1. - zmasq(:)
459
460	DO 99999 nsrf = 1, nbsrf
461
462	c chercher les indices:
463	DO j = 1, klon
464	ni(j) = 0
465	ENDDO
466	knon = 0
467	DO i = 1, klon
468
469	C pour determiner le domaine a traiter on utilise les surfaces "potentielles"
470	C
471	IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
472	knon = knon + 1
473	ni(knon) = i
474	ENDIF
475	ENDDO
476	c
477	if (check) THEN
478	write(,)'CLMAIN, nsrf, knon =',nsrf, knon
479	call flush(6)
480	endif
481	c
482	c variables pour avoir une sortie IOIPSL des INDEX
483	c
484	IF (debugindex) THEN
485	tabindx(:)=0.
486	c tabindx(1:knon)=(/FLOAT(i),i=1:knon/)
487	DO i=1,knon
488	tabindx(1:knon)=FLOAT(i)
489	END DO
490	debugtab(:,:)=0.
491	ndexbg(:)=0
492	CALL gath2cpl(tabindx,debugtab,klon,knon,iim,jjm,ni)
493	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1)
494	$ ,ndexbg)
495	ENDIF
496	IF (knon.EQ.0) GOTO 99999
497	DO j = 1, knon
498	i = ni(j)
499	ypct(j) = pctsrf(i,nsrf)
500	yts(j) = ts(i,nsrf)
501	cIM "slab" ocean
502	ytslab(i) = tslab(i)
503	c
504	ysnow(j) = snow(i,nsrf)
505	yqsurf(j) = qsurf(i,nsrf)
506	yalb(j) = albe(i,nsrf)
507	yalblw(j) = alblw(i,nsrf)
508	yrain_f(j) = rain_f(i)
509	ysnow_f(j) = snow_f(i)
510	yagesno(j) = agesno(i,nsrf)
511	yfder(j) = fder(i)
512	ytaux(j) = flux_u(i,1,nsrf)
513	ytauy(j) = flux_v(i,1,nsrf)
514	ysolsw(j) = solsw(i,nsrf)
515	ysollw(j) = sollw(i,nsrf)
516	ysollwdown(j) = sollwdown(i)
517	yrugos(j) = rugos(i,nsrf)
518	yrugoro(j) = rugoro(i)
519	yu1(j) = u1lay(i)
520	yv1(j) = v1lay(i)
521	yrads(j) = ysolsw(j)+ ysollw(j)
522	ypaprs(j,klev+1) = paprs(i,klev+1)
523	y_run_off_lic_0(j) = run_off_lic_0(i)
524	c -- LOOP
525	yu10mx(j) = u10m(i,nsrf)
526	yu10my(j) = v10m(i,nsrf)
527	ywindsp(j) = SQRT(yu10mx(j)yu10mx(j) + yu10my(j)yu10my(j) )
528	c -- LOOP
529	END DO
530	C
531	C IF bucket model for continent, copy soil water content
532	IF ( nsrf .eq. is_ter .and. .not. ok_veget ) THEN
533	DO j = 1, knon
534	i = ni(j)
535	yqsol(j) = qsol(i)
536	END DO
537	ELSE
538	yqsol(:)=0.
539	ENDIF
540	c$$$ PB ajour pour soil
541	DO k = 1, nsoilmx
542	DO j = 1, knon
543	i = ni(j)
544	ytsoil(j,k) = ftsoil(i,k,nsrf)
545	END DO
546	END DO
547	DO k = 1, klev
548	DO j = 1, knon
549	i = ni(j)
550	ypaprs(j,k) = paprs(i,k)
551	ypplay(j,k) = pplay(i,k)
552	ydelp(j,k) = delp(i,k)
553	yu(j,k) = u(i,k)
554	yv(j,k) = v(i,k)
555	yt(j,k) = t(i,k)
556	yq(j,k) = q(i,k)
557	ENDDO
558	ENDDO
559	c
560	c
561	c calculer Cdrag et les coefficients d'echange
562	CALL coefkz(nsrf, knon, ypaprs, ypplay,
563	cIM 261103
564	. ksta, ksta_ter,
565	cIM 261103
566	. yts, yrugos, yu, yv, yt, yq,
567	. yqsurf,
568	. ycoefm, ycoefh)
569	cIM 081204 BEG
570	cCR test
571	if (iflag_pbl.eq.1) then
572	cIM 081204 END
573	CALL coefkz2(nsrf, knon, ypaprs, ypplay,yt,
574	. ycoefm0, ycoefh0)
575	DO k = 1, klev
576	DO i = 1, knon
577	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
578	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
579	ENDDO
580	ENDDO
581	endif
582	c
583	cIM cf JLD : on seuille ycoefm et ycoefh
584	if (nsrf.eq.is_oce) then
585	do j=1,knon
586	c ycoefm(j,1)=min(ycoefm(j,1),1.1E-3)
587	ycoefm(j,1)=min(ycoefm(j,1),cdmmax)
588	c ycoefh(j,1)=min(ycoefh(j,1),1.1E-3)
589	ycoefh(j,1)=min(ycoefh(j,1),cdhmax)
590	enddo
591	endif
592
593	c
594	cIM: 261103
595	if (ok_kzmin) THEN
596	cIM cf FH: 201103 BEG
597	c Calcul d'une diffusion minimale pour les conditions tres stables.
598	call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm
599	. ,ycoefm0,ycoefh0)
600	c call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ ')
601	c call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN ')
602
603	if ( 1.eq.1 ) then
604	DO k = 1, klev
605	DO i = 1, knon
606	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
607	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
608	ENDDO
609	ENDDO
610	endif
611	cIM cf FH: 201103 END
612	endif !ok_kzmin
613	cIM: 261103
614
615
616	IF (iflag_pbl.ge.3) then
617
618	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
619	c MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
620	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
621
622	yzlay(1:knon,1)=
623	. RDyt(1:knon,1)/(0.5(ypaprs(1:knon,1)+ypplay(1:knon,1)))
624	. *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
625	do k=2,klev
626	yzlay(1:knon,k)=
627	. yzlay(1:knon,k-1)+RD0.5(yt(1:knon,k-1)+yt(1:knon,k))
628	. /ypaprs(1:knon,k)*(ypplay(1:knon,k-1)-ypplay(1:knon,k))/RG
629	enddo
630	do k=1,klev
631	yteta(1:knon,k)=
632	. yt(1:knon,k)(ypaprs(1:knon,1)/ypplay(1:knon,k))*rkappa
633	. (1.+0.61yq(1:knon,k))
634	enddo
635	yzlev(1:knon,1)=0.
636	yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
637	do k=2,klev
638	yzlev(1:knon,k)=0.5*(yzlay(1:knon,k)+yzlay(1:knon,k-1))
639	enddo
640	DO k = 1, klev+1
641	DO j = 1, knon
642	i = ni(j)
643	yq2(j,k)=q2(i,k,nsrf)
644	enddo
645	enddo
646
647
648	c Bug introduit volontairement pour converger avec les resultats
649	c du papier sur les thermiques.
650	if (1.eq.1) then
651	y_cd_m(1:knon) = ycoefm(1:knon,1)
652	y_cd_h(1:knon) = ycoefh(1:knon,1)
653	else
654	y_cd_h(1:knon) = ycoefm(1:knon,1)
655	y_cd_m(1:knon) = ycoefh(1:knon,1)
656	endif
657	call ustarhb(knon,yu,yv,y_cd_m, yustar)
658
659	if (prt_level > 9) THEN
660	WRITE(lunout,*)'USTAR = ',yustar
661	ENDIF
662
663	c iflag_pbl peut etre utilise comme longuer de melange
664
665	if (iflag_pbl.ge.11) then
666	call vdif_kcay(knon,dtime,rg,rd,ypaprs,yt
667	s ,yzlev,yzlay,yu,yv,yteta
668	s ,y_cd_m,yq2,q2diag,ykmm,ykmn,yustar,
669	s iflag_pbl)
670	else
671	call yamada4(knon,dtime,rg,rd,ypaprs,yt
672	s ,yzlev,yzlay,yu,yv,yteta
673	s ,y_cd_m,yq2,ykmm,ykmn,ykmq,yustar,
674	s iflag_pbl)
675	endif
676
677	ycoefm(1:knon,1)=y_cd_m(1:knon)
678	ycoefh(1:knon,1)=y_cd_h(1:knon)
679	ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
680	ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
681
682
683	ENDIF
684
685	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
686	c calculer la diffusion des vitesses "u" et "v"
687	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
688
689	CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
690	s y_d_u,y_flux_u)
691	CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
692	s y_d_v,y_flux_v)
693
694	c pour le couplage
695	ytaux = y_flux_u(:,1)
696	ytauy = y_flux_v(:,1)
697
698	c FH modif sur le cdrag temperature
699	c$$$PB : déplace dans clcdrag
700	c$$$ do i=1,knon
701	c$$$ ycoefh(i,1)=ycoefm(i,1)*0.8
702	c$$$ enddo
703
704	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
705	c calculer la diffusion de "q" et de "h"
706	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
707	CALL clqh(dtime, itap, date0,jour, debut,lafin,
708	e rlon, rlat, cufi, cvfi,
709	e knon, nsrf, ni, pctsrf,
710	e soil_model, ytsoil,yqsol,
711	e ok_veget, ocean, npas, nexca,
712	e rmu0, co2_ppm, yrugos, yrugoro,
713	e yu1, yv1, ycoefh,
714	e yt,yq,yts,ypaprs,ypplay,
715	e ydelp,yrads,yalb, yalblw, ysnow, yqsurf,
716	e yrain_f, ysnow_f, yfder, ytaux, ytauy,
717	c -- LOOP
718	e ywindsp,
719	c -- LOOP
720	c$$$ e ysollw, ysolsw,
721	e ysollw, ysollwdown, ysolsw,yfluxlat,
722	s pctsrf_new, yagesno,
723	s y_d_t, y_d_q, y_d_ts, yz0_new,
724	s y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,
725	s y_fqcalving,y_ffonte,y_run_off_lic_0,
726	cIM "slab" ocean
727	s y_flux_o, y_flux_g, ytslab, y_seaice)
728	c
729	c calculer la longueur de rugosite sur ocean
730	yrugm=0.
731	IF (nsrf.EQ.is_oce) THEN
732	DO j = 1, knon
733	yrugm(j) = 0.018ycoefm(j,1) (yu1(j)2+yv1(j)2)/RG
734	$ + 0.1114e-6 / sqrt(ycoefm(j,1) (yu1(j)2+yv1(j)2))
735	yrugm(j) = MAX(1.5e-05,yrugm(j))
736	ENDDO
737	ENDIF
738	DO j = 1, knon
739	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
740	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
741	yu1(j) = yu1(j) * ypct(j)
742	yv1(j) = yv1(j) * ypct(j)
743	ENDDO
744	c
745	DO k = 1, klev
746	DO j = 1, knon
747	i = ni(j)
748	ycoefh(j,k) = ycoefh(j,k) * ypct(j)
749	ycoefm(j,k) = ycoefm(j,k) * ypct(j)
750	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
751	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
752	C§§§ PB
753	flux_t(i,k,nsrf) = y_flux_t(j,k)
754	flux_q(i,k,nsrf) = y_flux_q(j,k)
755	flux_u(i,k,nsrf) = y_flux_u(j,k)
756	flux_v(i,k,nsrf) = y_flux_v(j,k)
757	c$$$ PB y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
758	c$$$ PB y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
759	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
760	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
761	c$$$ PB y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
762	c$$$ PB y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
763	ENDDO
764	ENDDO
765
766
767	evap(:,nsrf) = - flux_q(:,1,nsrf)
768	c
769	albe(:, nsrf) = 0.
770	alblw(:, nsrf) = 0.
771	snow(:, nsrf) = 0.
772	qsurf(:, nsrf) = 0.
773	rugos(:, nsrf) = 0.
774	fluxlat(:,nsrf) = 0.
775	DO j = 1, knon
776	i = ni(j)
777	d_ts(i,nsrf) = y_d_ts(j)
778	albe(i,nsrf) = yalb(j)
779	alblw(i,nsrf) = yalblw(j)
780	snow(i,nsrf) = ysnow(j)
781	qsurf(i,nsrf) = yqsurf(j)
782	rugos(i,nsrf) = yz0_new(j)
783	fluxlat(i,nsrf) = yfluxlat(j)
784	c$$$ pb rugmer(i) = yrugm(j)
785	IF (nsrf .EQ. is_oce) then
786	rugmer(i) = yrugm(j)
787	rugos(i,nsrf) = yrugm(j)
788	endif
789	cIM cf JLD ??
790	agesno(i,nsrf) = yagesno(j)
791	fqcalving(i,nsrf) = y_fqcalving(j)
792	ffonte(i,nsrf) = y_ffonte(j)
793	cdragh(i) = cdragh(i) + ycoefh(j,1)
794	cdragm(i) = cdragm(i) + ycoefm(j,1)
795	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
796	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
797	zu1(i) = zu1(i) + yu1(j)
798	zv1(i) = zv1(i) + yv1(j)
799	END DO
800	IF ( nsrf .eq. is_ter ) THEN
801	DO j = 1, knon
802	i = ni(j)
803	qsol(i) = yqsol(j)
804	END DO
805	END IF
806	IF ( nsrf .eq. is_lic ) THEN
807	DO j = 1, knon
808	i = ni(j)
809	run_off_lic_0(i) = y_run_off_lic_0(j)
810	END DO
811	END IF
812	c$$$ PB ajout pour soil
813	ftsoil(:,:,nsrf) = 0.
814	DO k = 1, nsoilmx
815	DO j = 1, knon
816	i = ni(j)
817	ftsoil(i, k, nsrf) = ytsoil(j,k)
818	END DO
819	END DO
820	c
821	#ifdef CRAY
822	DO k = 1, klev
823	DO j = 1, knon
824	i = ni(j)
825	#else
826	DO j = 1, knon
827	i = ni(j)
828	DO k = 1, klev
829	#endif
830	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
831	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
832	c$$$ PB flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
833	c$$$ flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
834	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
835	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
836	c$$$ PB flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
837	c$$$ flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
838	zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
839	ENDDO
840	ENDDO
841	c
842	c
843	#undef T2m
844	#define T2m
845	#ifdef T2m
846	ccc diagnostic t,q a 2m et u, v a 10m
847	c
848	DO j=1, knon
849	i = ni(j)
850	uzon(j) = yu(j,1) + y_d_u(j,1)
851	vmer(j) = yv(j,1) + y_d_v(j,1)
852	tair1(j) = yt(j,1) + y_d_t(j,1)
853	qair1(j) = yq(j,1) + y_d_q(j,1)
854	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1)))
855	& * (ypaprs(j,1)-ypplay(j,1))
856	tairsol(j) = yts(j) + y_d_ts(j)
857	rugo1(j) = yrugos(j)
858	IF(nsrf.EQ.is_oce) THEN
859	rugo1(j) = rugos(i,nsrf)
860	ENDIF
861	psfce(j)=ypaprs(j,1)
862	patm(j)=ypplay(j,1)
863	c
864	qairsol(j) = yqsurf(j)
865	ENDDO
866	c
867	CALL stdlevvar(klon, knon, nsrf, zxli,
868	& uzon, vmer, tair1, qair1, zgeo1,
869	& tairsol, qairsol, rugo1, psfce, patm,
870	cIM & yt2m, yq2m, yu10m)
871	& yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
872	cIM 081204 END
873	c
874	c
875	DO j=1, knon
876	i = ni(j)
877	t2m(i,nsrf)=yt2m(j)
878
879	c
880	q2m(i,nsrf)=yq2m(j)
881	c
882	c u10m, v10m : composantes du vent a 10m sans spirale de Ekman
883	u10m(i,nsrf)=(yu10m(j) * uzon(j))/sqrt(uzon(j)2+vmer(j)2)
884	v10m(i,nsrf)=(yu10m(j) * vmer(j))/sqrt(uzon(j)2+vmer(j)2)
885	c
886	ENDDO
887	c
888	cIM cf AM : pbl, HBTM
889	DO i = 1, knon
890	y_cd_h(i) = ycoefh(i,1)
891	y_cd_m(i) = ycoefm(i,1)
892	ENDDO
893	c print*,'appel hbtm2'
894	CALL HBTM(knon, ypaprs, ypplay,
895	. yt2m,yt10m,yq2m,yq10m,yustar,
896	. y_flux_t,y_flux_q,yu,yv,yt,yq,
897	. ypblh,ycapCL,yoliqCL,ycteiCL,ypblT,
898	. ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
899	c print*,'fin hbtm2'
900	c
901	DO j=1, knon
902	i = ni(j)
903	pblh(i,nsrf) = ypblh(j)
904	plcl(i,nsrf) = ylcl(j)
905	capCL(i,nsrf) = ycapCL(j)
906	oliqCL(i,nsrf) = yoliqCL(j)
907	cteiCL(i,nsrf) = ycteiCL(j)
908	pblT(i,nsrf) = ypblT(j)
909	therm(i,nsrf) = ytherm(j)
910	trmb1(i,nsrf) = ytrmb1(j)
911	trmb2(i,nsrf) = ytrmb2(j)
912	trmb3(i,nsrf) = ytrmb3(j)
913	ENDDO
914	c
915	#else
916	cIM cf AM : pbl, HBTM
917	DO j=1, knon
918	i = ni(j)
919	pblh(i,nsrf) = 0.
920	plcl(i,nsrf) = 0.
921	capCL(i,nsrf) = 0.
922	oliqCL(i,nsrf) = 0.
923	cteiCL(i,nsrf) = 0.
924	pblT(i,nsrf) = 0.
925	therm(i,nsrf) = 0.
926	trmb1(i,nsrf) = 0.
927	trmb2(i,nsrf) = 0.
928	trmb3(i,nsrf) = 0.
929	ENDDO
930	DO j=1, knon
931	i = ni(j)
932	t2m(i,nsrf)=0.
933	q2m(i,nsrf)=0.
934	u10m(i,nsrf)=0.
935	v10m(i,nsrf)=0.
936	ENDDO
937	#endif
938
939	do j=1,knon
940	do k=1,klev+1
941	i=ni(j)
942	q2(i,k,nsrf)=yq2(j,k)
943	enddo
944	enddo
945	cIM "slab" ocean
946	IF(OCEAN.EQ.'slab '.OR.OCEAN.EQ.'force ') THEN
947	IF (nsrf.EQ.is_oce) THEN
948	DO j = 1, knon
949	c on projette sur la grille globale
950	i = ni(j)
951	IF(pctsrf_new(i,is_oce).GT.epsfra) THEN
952	flux_o(i) = y_flux_o(j)
953	ELSE
954	flux_o(i) = 0.
955	ENDIF
956	ENDDO
957	ENDIF
958	c
959	IF (nsrf.EQ.is_sic) THEN
960	DO j = 1, knon
961	i = ni(j)
962	cIM 230604 on pondere lorsque l'on fait le bilan au sol : flux_g(i) = y_flux_g(j)*ypct(j)
963	IF(pctsrf_new(i,is_sic).GT.epsfra) THEN
964	flux_g(i) = y_flux_g(j)
965	ELSE
966	flux_g(i) = 0.
967	ENDIF
968	ENDDO
969	ENDIF !nsrf.EQ.is_sic
970	ENDIF !OCEAN
971	c
972	IF(OCEAN.EQ.'slab ') THEN
973	IF(nsrf.EQ.is_oce) then
974	tslab(1:klon) = ytslab(1:klon)
975	seaice(1:klon) = y_seaice(1:klon)
976	ENDIF !nsrf
977	ENDIF !OCEAN
978	99999 CONTINUE
979	C
980	C On utilise les nouvelles surfaces
981	C A rajouter: conservation de l'albedo
982	C
983	rugos(:,is_oce) = rugmer
984	pctsrf = pctsrf_new
985
986	RETURN
987	END
988	SUBROUTINE clqh(dtime,itime, date0,jour,debut,lafin,
989	e rlon, rlat, cufi, cvfi,
990	e knon, nisurf, knindex, pctsrf,
991	$ soil_model,tsoil,qsol,
992	e ok_veget, ocean, npas, nexca,
993	e rmu0, co2_ppm, rugos, rugoro,
994	e u1lay,v1lay,coef,
995	e t,q,ts,paprs,pplay,
996	e delp,radsol,albedo,alblw,snow,qsurf,
997	e precip_rain, precip_snow, fder, taux, tauy,
998	c -- LOOP
999	e ywindsp,
1000	c -- LOOP
1001	$ sollw, sollwdown, swnet,fluxlat,
1002	s pctsrf_new, agesno,
1003	s d_t, d_q, d_ts, z0_new,
1004	s flux_t, flux_q,dflux_s,dflux_l,
1005	s fqcalving,ffonte,run_off_lic_0,
1006	cIM "slab" ocean
1007	s flux_o,flux_g,tslab,seaice)
1008
1009	USE interface_surf
1010
1011	IMPLICIT none
1012	c======================================================================
1013	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
1014	c Objet: diffusion verticale de "q" et de "h"
1015	c======================================================================
1016	#include "dimensions.h"
1017	#include "dimphy.h"
1018	#include "YOMCST.h"
1019	#include "YOETHF.h"
1020	#include "FCTTRE.h"
1021	#include "indicesol.h"
1022	#include "dimsoil.h"
1023	#include "iniprint.h"
1024
1025	c Arguments:
1026	INTEGER knon
1027	REAL dtime ! intervalle du temps (s)
1028	real date0
1029	REAL u1lay(klon) ! vitesse u de la 1ere couche (m/s)
1030	REAL v1lay(klon) ! vitesse v de la 1ere couche (m/s)
1031	REAL coef(klon,klev) ! le coefficient d'echange (m**2/s)
1032	c multiplie par le cisaillement du
1033	c vent (dV/dz); la premiere valeur
1034	c indique la valeur de Cdrag (sans unite)
1035	REAL t(klon,klev) ! temperature (K)
1036	REAL q(klon,klev) ! humidite specifique (kg/kg)
1037	REAL ts(klon) ! temperature du sol (K)
1038	REAL evap(klon) ! evaporation au sol
1039	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
1040	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
1041	REAL delp(klon,klev) ! epaisseur de couche en pression (Pa)
1042	REAL radsol(klon) ! ray. net au sol (Solaire+IR) W/m2
1043	REAL albedo(klon) ! albedo de la surface
1044	REAL alblw(klon)
1045	REAL snow(klon) ! hauteur de neige
1046	REAL qsurf(klon) ! humidite de l'air au dessus de la surface
1047	real precip_rain(klon), precip_snow(klon)
1048	REAL agesno(klon)
1049	REAL rugoro(klon)
1050	REAL run_off_lic_0(klon)! runof glacier au pas de temps precedent
1051	integer jour ! jour de l'annee en cours
1052	real rmu0(klon) ! cosinus de l'angle solaire zenithal
1053	real rugos(klon) ! rugosite
1054	integer knindex(klon)
1055	real pctsrf(klon,nbsrf)
1056	real rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
1057	logical ok_veget
1058	REAL co2_ppm ! taux CO2 atmosphere
1059	character*6 ocean
1060	integer npas, nexca
1061	c -- LOOP
1062	REAL yu10mx(klon)
1063	REAL yu10my(klon)
1064	REAL ywindsp(klon)
1065	c -- LOOP
1066
1067
1068	c
1069	REAL d_t(klon,klev) ! incrementation de "t"
1070	REAL d_q(klon,klev) ! incrementation de "q"
1071	REAL d_ts(klon) ! incrementation de "ts"
1072	REAL flux_t(klon,klev) ! (diagnostic) flux de la chaleur
1073	c sensible, flux de Cp*T, positif vers
1074	c le bas: j/(m**2 s) c.a.d.: W/m2
1075	REAL flux_q(klon,klev) ! flux de la vapeur d'eau:kg/(m**2 s)
1076	REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
1077	REAL dflux_l(klon) ! derivee du flux latent dF/dTs
1078	cIM cf JLD
1079	c Flux thermique utiliser pour fondre la neige
1080	REAL ffonte(klon)
1081	c Flux d'eau "perdue" par la surface et nécessaire pour que limiter la
1082	c hauteur de neige, en kg/m2/s
1083	REAL fqcalving(klon)
1084	cIM "slab" ocean
1085	REAL tslab(klon) !temperature du slab ocean (K) (OCEAN='slab ')
1086	REAL seaice(klon) ! glace de mer en kg/m2
1087	REAL flux_o(klon) ! flux entre l'ocean et l'atmosphere W/m2
1088	REAL flux_g(klon) ! flux entre l'ocean et la glace de mer W/m2
1089	c
1090	c======================================================================
1091	REAL t_grnd ! temperature de rappel pour glace de mer
1092	PARAMETER (t_grnd=271.35)
1093	REAL t_coup
1094	PARAMETER(t_coup=273.15)
1095	c======================================================================
1096	INTEGER i, k
1097	REAL zx_cq(klon,klev)
1098	REAL zx_dq(klon,klev)
1099	REAL zx_ch(klon,klev)
1100	REAL zx_dh(klon,klev)
1101	REAL zx_buf1(klon)
1102	REAL zx_buf2(klon)
1103	REAL zx_coef(klon,klev)
1104	REAL local_h(klon,klev) ! enthalpie potentielle
1105	REAL local_q(klon,klev)
1106	REAL local_ts(klon)
1107	REAL psref(klon) ! pression de reference pour temperature potent.
1108	REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
1109	c======================================================================
1110	c contre-gradient pour la vapeur d'eau: (kg/kg)/metre
1111	REAL gamq(klon,2:klev)
1112	c contre-gradient pour la chaleur sensible: Kelvin/metre
1113	REAL gamt(klon,2:klev)
1114	REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)
1115	REAL zdelz
1116	c======================================================================
1117	#include "compbl.h"
1118	c======================================================================
1119	c Rajout pour l'interface
1120	integer itime
1121	integer nisurf
1122	logical debut, lafin
1123	real zlev1(klon)
1124	real fder(klon), taux(klon), tauy(klon)
1125	real temp_air(klon), spechum(klon)
1126	real epot_air(klon), ccanopy(klon)
1127	real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
1128	real petBcoef(klon), peqBcoef(klon)
1129	real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
1130	real p1lay(klon)
1131	c$$$C PB ajout pour soil
1132	LOGICAL soil_model
1133	REAL tsoil(klon, nsoilmx)
1134	REAL qsol(klon)
1135
1136	! Parametres de sortie
1137	real fluxsens(klon), fluxlat(klon)
1138	real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
1139	real emis_new(klon), z0_new(klon)
1140	real pctsrf_new(klon,nbsrf)
1141	c JLD
1142	real zzpk
1143	C
1144	character (len = 20) :: modname = 'Debut clqh'
1145	LOGICAL check
1146	PARAMETER (check=.false.)
1147	C
1148	if (check) THEN
1149	write(,) modname,' nisurf=',nisurf
1150	call flush(6)
1151	endif
1152	c
1153	if (check) THEN
1154	WRITE(,)' qsurf (min, max)'
1155	$ , minval(qsurf(1:knon)), maxval(qsurf(1:knon))
1156	call flush(6)
1157	ENDIF
1158	C
1159	C
1160	if (iflag_pbl.eq.1) then
1161	do k = 3, klev
1162	do i = 1, knon
1163	gamq(i,k)= 0.0
1164	gamt(i,k)= -1.0e-03
1165	enddo
1166	enddo
1167	do i = 1, knon
1168	gamq(i,2) = 0.0
1169	gamt(i,2) = -2.5e-03
1170	enddo
1171	else
1172	do k = 2, klev
1173	do i = 1, knon
1174	gamq(i,k) = 0.0
1175	gamt(i,k) = 0.0
1176	enddo
1177	enddo
1178	endif
1179
1180	DO i = 1, knon
1181	psref(i) = paprs(i,1) !pression de reference est celle au sol
1182	local_ts(i) = ts(i)
1183	ENDDO
1184	DO k = 1, klev
1185	DO i = 1, knon
1186	zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
1187	zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
1188	local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)
1189	local_q(i,k) = q(i,k)
1190	ENDDO
1191	ENDDO
1192	c
1193	c Convertir les coefficients en variables convenables au calcul:
1194	c
1195	c
1196	DO k = 2, klev
1197	DO i = 1, knon
1198	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
1199	. (paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
1200	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
1201	ENDDO
1202	ENDDO
1203	c
1204	c Preparer les flux lies aux contre-gardients
1205	c
1206	DO k = 2, klev
1207	DO i = 1, knon
1208	zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k)
1209	. *(pplay(i,k-1)-pplay(i,k))
1210	z_gamaq(i,k) = gamq(i,k) * zdelz
1211	z_gamah(i,k) = gamt(i,k) * zdelz RCPD zx_pkh(i,k)
1212	ENDDO
1213	ENDDO
1214	DO i = 1, knon
1215	zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
1216	zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev)
1217	. -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)
1218	zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
1219	c
1220	zzpk=(pplay(i,klev)/psref(i))**RKAPPA
1221	zx_buf2(i) = zzpk*delp(i,klev) + zx_coef(i,klev)
1222	zx_ch(i,klev) = (local_h(i,klev)zzpkdelp(i,klev)
1223	. -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
1224	zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
1225	ENDDO
1226	DO k = klev-1, 2 , -1
1227	DO i = 1, knon
1228	zx_buf1(i) = delp(i,k)+zx_coef(i,k)
1229	. +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))
1230	zx_cq(i,k) = (local_q(i,k)*delp(i,k)
1231	. +zx_coef(i,k+1)*zx_cq(i,k+1)
1232	. +zx_coef(i,k+1)*z_gamaq(i,k+1)
1233	. -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)
1234	zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
1235	c
1236	zzpk=(pplay(i,k)/psref(i))**RKAPPA
1237	zx_buf2(i) = zzpk*delp(i,k)+zx_coef(i,k)
1238	. +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
1239	zx_ch(i,k) = (local_h(i,k)zzpkdelp(i,k)
1240	. +zx_coef(i,k+1)*zx_ch(i,k+1)
1241	. +zx_coef(i,k+1)*z_gamah(i,k+1)
1242	. -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
1243	zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
1244	ENDDO
1245	ENDDO
1246	C
1247	C nouvelle formulation JL Dufresne
1248	C
1249	C q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt
1250	C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
1251	C
1252	DO i = 1, knon
1253	zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))
1254	zx_cq(i,1) = (local_q(i,1)*delp(i,1)
1255	. +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2)))
1256	. /zx_buf1(i)
1257	zx_dq(i,1) = -1. * RG / zx_buf1(i)
1258	c
1259	zzpk=(pplay(i,1)/psref(i))**RKAPPA
1260	zx_buf2(i) = zzpkdelp(i,1) + zx_coef(i,2)(1.-zx_dh(i,2))
1261	zx_ch(i,1) = (local_h(i,1)zzpkdelp(i,1)
1262	. +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
1263	. /zx_buf2(i)
1264	zx_dh(i,1) = -1. * RG / zx_buf2(i)
1265	ENDDO
1266
1267	C Appel a interfsurf (appel generique) routine d'interface avec la surface
1268
1269	c initialisation
1270	petAcoef =0.
1271	peqAcoef = 0.
1272	petBcoef =0.
1273	peqBcoef = 0.
1274	p1lay =0.
1275
1276	c do i = 1, knon
1277	petAcoef(1:knon) = zx_ch(1:knon,1)
1278	peqAcoef(1:knon) = zx_cq(1:knon,1)
1279	petBcoef(1:knon) = zx_dh(1:knon,1)
1280	peqBcoef(1:knon) = zx_dq(1:knon,1)
1281	tq_cdrag(1:knon) =coef(1:knon,1)
1282	temp_air(1:knon) =t(1:knon,1)
1283	epot_air(1:knon) =local_h(1:knon,1)
1284	spechum(1:knon)=q(1:knon,1)
1285	p1lay(1:knon) = pplay(1:knon,1)
1286	zlev1(1:knon) = delp(1:knon,1)
1287	c swnet = swdown * (1. - albedo)
1288	c
1289	cIM swdown=flux SW incident sur terres
1290	cIM swdown=flux SW net sur les autres surfaces
1291	cIM swdown(1:knon) = swnet(1:knon)
1292	if(nisurf.eq.is_ter) THEN
1293	swdown(1:knon) = swnet(1:knon)/(1-albedo(1:knon))
1294	else
1295	swdown(1:knon) = swnet(1:knon)
1296	endif
1297	c enddo
1298	ccanopy = co2_ppm
1299
1300	CALL interfsurf(itime, dtime, date0, jour, rmu0,
1301	e klon, iim, jjm, nisurf, knon, knindex, pctsrf,
1302	e rlon, rlat, cufi, cvfi,
1303	e debut, lafin, ok_veget, soil_model, nsoilmx,tsoil, qsol,
1304	e zlev1, u1lay, v1lay, temp_air, spechum, epot_air, ccanopy,
1305	e tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
1306	e precip_rain, precip_snow, sollw, sollwdown, swnet, swdown,
1307	e fder, taux, tauy,
1308	c -- LOOP
1309	e ywindsp,
1310	c -- LOOP
1311	e rugos, rugoro,
1312	e albedo, snow, qsurf,
1313	e ts, p1lay, psref, radsol,
1314	e ocean, npas, nexca, zmasq,
1315	s evap, fluxsens, fluxlat, dflux_l, dflux_s,
1316	s tsol_rad, tsurf_new, alb_new, alblw, emis_new, z0_new,
1317	s pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0,
1318	cIM "slab" ocean
1319	s flux_o, flux_g, tslab, seaice)
1320
1321
1322	do i = 1, knon
1323	flux_t(i,1) = fluxsens(i)
1324	flux_q(i,1) = - evap(i)
1325	d_ts(i) = tsurf_new(i) - ts(i)
1326	albedo(i) = alb_new(i)
1327	enddo
1328
1329	c==== une fois on a zx_h_ts, on peut faire l'iteration ========
1330	DO i = 1, knon
1331	local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)flux_t(i,1)dtime
1332	local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)flux_q(i,1)dtime
1333	ENDDO
1334	DO k = 2, klev
1335	DO i = 1, knon
1336	local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)
1337	local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
1338	ENDDO
1339	ENDDO
1340	c======================================================================
1341	c== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
1342	c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
1343	DO k = 2, klev
1344	DO i = 1, knon
1345	flux_q(i,k) = (zx_coef(i,k)/RG/dtime)
1346	. * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))
1347	flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
1348	. * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
1349	. / zx_pkh(i,k)
1350	ENDDO
1351	ENDDO
1352	c======================================================================
1353	C Calcul tendances
1354	DO k = 1, klev
1355	DO i = 1, knon
1356	d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)
1357	d_q(i,k) = local_q(i,k) - q(i,k)
1358	ENDDO
1359	ENDDO
1360	c
1361
1362	RETURN
1363	END
1364	SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
1365	e paprs,pplay,delp,
1366	s d_ven,flux_v)
1367	IMPLICIT none
1368	c======================================================================
1369	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
1370	c Objet: diffusion vertical de la vitesse "ven"
1371	c======================================================================
1372	c Arguments:
1373	c dtime----input-R- intervalle du temps (en second)
1374	c u1lay----input-R- vent u de la premiere couche (m/s)
1375	c v1lay----input-R- vent v de la premiere couche (m/s)
1376	c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
1377	c le cisaillement du vent (dV/dz); la premiere
1378	c valeur indique la valeur de Cdrag (sans unite)
1379	c t--------input-R- temperature (K)
1380	c ven------input-R- vitesse horizontale (m/s)
1381	c paprs----input-R- pression a inter-couche (Pa)
1382	c pplay----input-R- pression au milieu de couche (Pa)
1383	c delp-----input-R- epaisseur de couche (Pa)
1384	c
1385	c
1386	c d_ven----output-R- le changement de "ven"
1387	c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
1388	c======================================================================
1389	#include "dimensions.h"
1390	#include "dimphy.h"
1391	#include "iniprint.h"
1392	INTEGER knon
1393	REAL dtime
1394	REAL u1lay(klon), v1lay(klon)
1395	REAL coef(klon,klev)
1396	REAL t(klon,klev), ven(klon,klev)
1397	REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
1398	REAL d_ven(klon,klev)
1399	REAL flux_v(klon,klev)
1400	c======================================================================
1401	#include "YOMCST.h"
1402	c======================================================================
1403	INTEGER i, k
1404	REAL zx_cv(klon,2:klev)
1405	REAL zx_dv(klon,2:klev)
1406	REAL zx_buf(klon)
1407	REAL zx_coef(klon,klev)
1408	REAL local_ven(klon,klev)
1409	REAL zx_alf1(klon), zx_alf2(klon)
1410	c======================================================================
1411	DO k = 1, klev
1412	DO i = 1, knon
1413	local_ven(i,k) = ven(i,k)
1414	ENDDO
1415	ENDDO
1416	c======================================================================
1417	DO i = 1, knon
1418	ccc zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
1419	zx_alf1(i) = 1.0
1420	zx_alf2(i) = 1.0 - zx_alf1(i)
1421	zx_coef(i,1) = coef(i,1)
1422	. * (1.0+SQRT(u1lay(i)2+v1lay(i)2))
1423	. * pplay(i,1)/(RD*t(i,1))
1424	zx_coef(i,1) = zx_coef(i,1) * dtime*RG
1425	ENDDO
1426	c======================================================================
1427	DO k = 2, klev
1428	DO i = 1, knon
1429	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
1430	. (paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
1431	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
1432	ENDDO
1433	ENDDO
1434	c======================================================================
1435	DO i = 1, knon
1436	zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
1437	zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
1438	zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
1439	. /zx_buf(i)
1440	ENDDO
1441	DO k = 3, klev
1442	DO i = 1, knon
1443	zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
1444	. + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
1445	zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
1446	. +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
1447	zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
1448	ENDDO
1449	ENDDO
1450	DO i = 1, knon
1451	local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
1452	. +zx_coef(i,klev)*zx_cv(i,klev) )
1453	. / ( delp(i,klev) + zx_coef(i,klev)
1454	. -zx_coef(i,klev)*zx_dv(i,klev) )
1455	ENDDO
1456	DO k = klev-1, 1, -1
1457	DO i = 1, knon
1458	local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
1459	ENDDO
1460	ENDDO
1461	c======================================================================
1462	c== flux_v est le flux de moment angulaire (positif vers bas)
1463	c== dont l'unite est: (kg m/s)/(m**2 s)
1464	DO i = 1, knon
1465	flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
1466	. (local_ven(i,1)zx_alf1(i)
1467	. +local_ven(i,2)*zx_alf2(i))
1468	ENDDO
1469	DO k = 2, klev
1470	DO i = 1, knon
1471	flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
1472	. * (local_ven(i,k)-local_ven(i,k-1))
1473	ENDDO
1474	ENDDO
1475	c
1476	DO k = 1, klev
1477	DO i = 1, knon
1478	d_ven(i,k) = local_ven(i,k) - ven(i,k)
1479	ENDDO
1480	ENDDO
1481	c
1482	RETURN
1483	END
1484	SUBROUTINE coefkz(nsrf, knon, paprs, pplay,
1485	cIM 261103
1486	. ksta, ksta_ter,
1487	cIM 261103
1488	. ts, rugos,
1489	. u,v,t,q,
1490	. qsurf,
1491	. pcfm, pcfh)
1492	IMPLICIT none
1493	c======================================================================
1494	c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
1495	c (une version strictement identique a l'ancien modele)
1496	c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
1497	c coefficients d'echange turbulent dans l'atmosphere.
1498	c Arguments:
1499	c nsrf-----input-I- indicateur de la nature du sol
1500	c knon-----input-I- nombre de points a traiter
1501	c paprs----input-R- pression a chaque intercouche (en Pa)
1502	c pplay----input-R- pression au milieu de chaque couche (en Pa)
1503	c ts-------input-R- temperature du sol (en Kelvin)
1504	c rugos----input-R- longeur de rugosite (en m)
1505	c u--------input-R- vitesse u
1506	c v--------input-R- vitesse v
1507	c t--------input-R- temperature (K)
1508	c q--------input-R- vapeur d'eau (kg/kg)
1509	c
1510	c itop-----output-I- numero de couche du sommet de la couche limite
1511	c pcfm-----output-R- coefficients a calculer (vitesse)
1512	c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1513	c======================================================================
1514	#include "dimensions.h"
1515	#include "dimphy.h"
1516	#include "YOMCST.h"
1517	#include "indicesol.h"
1518	#include "iniprint.h"
1519	c
1520	c Arguments:
1521	c
1522	INTEGER knon, nsrf
1523	REAL ts(klon)
1524	REAL paprs(klon,klev+1), pplay(klon,klev)
1525	REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
1526	REAL rugos(klon)
1527	c
1528	REAL pcfm(klon,klev), pcfh(klon,klev)
1529	INTEGER itop(klon)
1530	c
1531	c Quelques constantes et options:
1532	c
1533	REAL cepdu2, ckap, cb, cc, cd, clam
1534	PARAMETER (cepdu2 =(0.1)**2)
1535	PARAMETER (CKAP=0.4)
1536	PARAMETER (cb=5.0)
1537	PARAMETER (cc=5.0)
1538	PARAMETER (cd=5.0)
1539	PARAMETER (clam=160.0)
1540	REAL ratqs ! largeur de distribution de vapeur d'eau
1541	PARAMETER (ratqs=0.05)
1542	LOGICAL richum ! utilise le nombre de Richardson humide
1543	PARAMETER (richum=.TRUE.)
1544	REAL ric ! nombre de Richardson critique
1545	PARAMETER(ric=0.4)
1546	REAL prandtl
1547	PARAMETER (prandtl=0.4)
1548	REAL kstable ! diffusion minimale (situation stable)
1549	! GKtest
1550	! PARAMETER (kstable=1.0e-10)
1551	REAL ksta, ksta_ter
1552	cIM: 261103 REAL kstable_ter, kstable_sinon
1553	cIM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6)
1554	cIM: 261103 PARAMETER (kstable_ter = 1.0e-8)
1555	cIM: 261103 PARAMETER (kstable_ter = 1.0e-10)
1556	cIM: 261103 PARAMETER (kstable_sinon = 1.0e-10)
1557	! fin GKtest
1558	REAL mixlen ! constante controlant longueur de melange
1559	PARAMETER (mixlen=35.0)
1560	INTEGER isommet ! le sommet de la couche limite
1561	PARAMETER (isommet=klev)
1562	LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
1563	PARAMETER (tvirtu=.TRUE.)
1564	LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
1565	PARAMETER (opt_ec=.FALSE.)
1566
1567	#include "compbl.h"
1568	c
1569	c Variables locales:
1570	c
1571	INTEGER i, k, kk !IM 120704
1572	REAL zgeop(klon,klev)
1573	REAL zmgeom(klon)
1574	REAL zri(klon)
1575	REAL zl2(klon)
1576
1577	REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
1578	REAL pcfm1(klon), pcfh1(klon)
1579	c
1580	REAL zdphi, zdu2, ztvd, ztvu, zcdn
1581	REAL zscf
1582	REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
1583	REAL z2geomf, zalh2, zalm2, zscfh, zscfm
1584	REAL t_coup
1585	PARAMETER (t_coup=273.15)
1586	cIM
1587	LOGICAL check
1588	PARAMETER (check=.false.)
1589	c
1590	c contre-gradient pour la chaleur sensible: Kelvin/metre
1591	REAL gamt(2:klev)
1592	real qsurf(klon)
1593	c
1594	LOGICAL appel1er
1595	SAVE appel1er
1596	c
1597	c Fonctions thermodynamiques et fonctions d'instabilite
1598	REAL fsta, fins, x
1599	LOGICAL zxli ! utiliser un jeu de fonctions simples
1600	PARAMETER (zxli=.FALSE.)
1601	c
1602	#include "YOETHF.h"
1603	#include "FCTTRE.h"
1604	fsta(x) = 1.0 / (1.0+10.0x(1+8.0*x))
1605	fins(x) = SQRT(1.0-18.0*x)
1606	c
1607	DATA appel1er /.TRUE./
1608	c
1609	IF (appel1er) THEN
1610	if (prt_level > 9) THEN
1611	WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
1612	WRITE(lunout,*)'coefkz, richum:', richum
1613	IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
1614	WRITE(lunout,*)'coefkz, isommet:', isommet
1615	WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
1616	appel1er = .FALSE.
1617	endif
1618	ENDIF
1619	c
1620	c Initialiser les sorties
1621	c
1622	DO k = 1, klev
1623	DO i = 1, knon
1624	pcfm(i,k) = 0.0
1625	pcfh(i,k) = 0.0
1626	ENDDO
1627	ENDDO
1628	DO i = 1, knon
1629	itop(i) = 0
1630	ENDDO
1631
1632	c
1633	c Prescrire la valeur de contre-gradient
1634	c
1635	if (iflag_pbl.eq.1) then
1636	DO k = 3, klev
1637	gamt(k) = -1.0E-03
1638	ENDDO
1639	gamt(2) = -2.5E-03
1640	else
1641	DO k = 2, klev
1642	gamt(k) = 0.0
1643	ENDDO
1644	ENDIF
1645	cIM cf JLD/ GKtest
1646	IF ( nsrf .NE. is_oce ) THEN
1647	cIM 261103 kstable = kstable_ter
1648	kstable = ksta_ter
1649	ELSE
1650	cIM 261103 kstable = kstable_sinon
1651	kstable = ksta
1652	ENDIF
1653	cIM cf JLD/ GKtest fin
1654	c
1655	c Calculer les geopotentiels de chaque couche
1656	c
1657	DO i = 1, knon
1658	zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1659	. * (paprs(i,1)-pplay(i,1))
1660	ENDDO
1661	DO k = 2, klev
1662	DO i = 1, knon
1663	zgeop(i,k) = zgeop(i,k-1)
1664	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1665	. * (pplay(i,k-1)-pplay(i,k))
1666	ENDDO
1667	ENDDO
1668	c
1669	c Calculer le frottement au sol (Cdrag)
1670	c
1671	DO i = 1, knon
1672	u1(i) = u(i,1)
1673	v1(i) = v(i,1)
1674	t1(i) = t(i,1)
1675	q1(i) = q(i,1)
1676	z1(i) = zgeop(i,1)
1677	ENDDO
1678	c
1679	CALL clcdrag(klon, knon, nsrf, zxli,
1680	$ u1, v1, t1, q1, z1,
1681	$ ts, qsurf, rugos,
1682	$ pcfm1, pcfh1)
1683	cIM $ ts, qsurf, rugos,
1684	C
1685	DO i = 1, knon
1686	pcfm(i,1)=pcfm1(i)
1687	pcfh(i,1)=pcfh1(i)
1688	ENDDO
1689	c
1690	c Calculer les coefficients turbulents dans l'atmosphere
1691	c
1692	DO i = 1, knon
1693	itop(i) = isommet
1694	ENDDO
1695
1696
1697	DO k = 2, isommet
1698	DO i = 1, knon
1699	zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
1700	. +(v(i,k)-v(i,k-1))**2)
1701	zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
1702	zdphi =zmgeom(i) / 2.0
1703	zt = (t(i,k)+t(i,k-1)) * 0.5
1704	zq = (q(i,k)+q(i,k-1)) * 0.5
1705	c
1706	c calculer Qs et dQs/dT:
1707	c
1708	IF (thermcep) THEN
1709	zdelta = MAX(0.,SIGN(1.,RTT-zt))
1710	zcvm5 = R5LESRLVTT/RCPD/(1.0+RVTMP2zq)*(1.-zdelta)
1711	. + R5IESRLSTT/RCPD/(1.0+RVTMP2zq)*zdelta
1712	zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
1713	zqs = MIN(0.5,zqs)
1714	zcor = 1./(1.-RETV*zqs)
1715	zqs = zqs*zcor
1716	zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
1717	ELSE
1718	IF (zt .LT. t_coup) THEN
1719	zqs = qsats(zt) / pplay(i,k)
1720	zdqs = dqsats(zt,zqs)
1721	ELSE
1722	zqs = qsatl(zt) / pplay(i,k)
1723	zdqs = dqsatl(zt,zqs)
1724	ENDIF
1725	ENDIF
1726	c
1727	c calculer la fraction nuageuse (processus humide):
1728	c
1729	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
1730	zfr = MAX(0.0,MIN(1.0,zfr))
1731	IF (.NOT.richum) zfr = 0.0
1732	c
1733	c calculer le nombre de Richardson:
1734	c
1735	IF (tvirtu) THEN
1736	ztvd =( t(i,k)
1737	. + zdphi/RCPD/(1.+RVTMP2*zq)
1738	. ( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
1739	. )(1.+RETVq(i,k))
1740	ztvu =( t(i,k-1)
1741	. - zdphi/RCPD/(1.+RVTMP2*zq)
1742	. ( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
1743	. )(1.+RETVq(i,k-1))
1744	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
1745	zri(i) = zri(i)
1746	. + zmgeom(i)zmgeom(i)/RGgamt(k)
1747	. (paprs(i,k)/101325.0)*RKAPPA
1748	. /(zdu20.5(ztvd+ztvu))
1749	c
1750	ELSE ! calcul de Ridchardson compatible LMD5
1751	c
1752	zri(i) =(RCPD*(t(i,k)-t(i,k-1))
1753	. -RD0.5(t(i,k)+t(i,k-1))/paprs(i,k)
1754	. *(pplay(i,k)-pplay(i,k-1))
1755	. )zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
1756	zri(i) = zri(i) +
1757	. zmgeom(i)zmgeom(i)gamt(k)/RG
1758	cSB . /(paprs(i,k)/101325.0)**RKAPPA
1759	. (paprs(i,k)/101325.0)*RKAPPA
1760	. /(zdu20.5(t(i,k-1)+t(i,k)))
1761	ENDIF
1762	c
1763	c finalement, les coefficients d'echange sont obtenus:
1764	c
1765	zcdn=SQRT(zdu2) / zmgeom(i) * RG
1766	c
1767	IF (opt_ec) THEN
1768	z2geomf=zgeop(i,k-1)+zgeop(i,k)
1769	zalm2=(0.5ckap/RGz2geomf
1770	. /(1.+0.5ckap/rg/clamz2geomf))**2
1771	zalh2=(0.5ckap/rgz2geomf
1772	. /(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
1773	IF (zri(i).LT.0.0) THEN ! situation instable
1774	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3
1775	. / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
1776	zscf = SQRT(-zri(i)*zscf)
1777	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
1778	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
1779	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
1780	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
1781	ELSE ! situation stable
1782	zscf=SQRT(1.+cd*zri(i))
1783	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
1784	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
1785	ENDIF
1786	ELSE
1787	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
1788	. /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
1789	pcfm(i,k)=sqrt(max(zcdnzcdn(ric-zri(i))/ric, kstable))
1790	pcfm(i,k)= zl2(i)* pcfm(i,k)
1791	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1792	ENDIF
1793	ENDDO
1794	ENDDO
1795	c
1796	c Au-dela du sommet, pas de diffusion turbulente:
1797	c
1798	DO i = 1, knon
1799	IF (itop(i)+1 .LE. klev) THEN
1800	DO k = itop(i)+1, klev
1801	pcfh(i,k) = 0.0
1802	pcfm(i,k) = 0.0
1803	ENDDO
1804	ENDIF
1805	ENDDO
1806	c
1807	RETURN
1808	END
1809
1810	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t,
1811	. pcfm, pcfh)
1812	IMPLICIT none
1813	c======================================================================
1814	c J'introduit un peu de diffusion sauf dans les endroits
1815	c ou une forte inversion est presente
1816	c On peut dire qu'il represente la convection peu profonde
1817	c
1818	c Arguments:
1819	c nsrf-----input-I- indicateur de la nature du sol
1820	c knon-----input-I- nombre de points a traiter
1821	c paprs----input-R- pression a chaque intercouche (en Pa)
1822	c pplay----input-R- pression au milieu de chaque couche (en Pa)
1823	c t--------input-R- temperature (K)
1824	c
1825	c pcfm-----output-R- coefficients a calculer (vitesse)
1826	c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1827	c======================================================================
1828	#include "dimensions.h"
1829	#include "dimphy.h"
1830	#include "YOMCST.h"
1831	#include "indicesol.h"
1832	#include "iniprint.h"
1833	c
1834	c Arguments:
1835	c
1836	INTEGER knon, nsrf
1837	REAL paprs(klon,klev+1), pplay(klon,klev)
1838	REAL t(klon,klev)
1839	c
1840	REAL pcfm(klon,klev), pcfh(klon,klev)
1841	c
1842	c Quelques constantes et options:
1843	c
1844	REAL prandtl
1845	PARAMETER (prandtl=0.4)
1846	REAL kstable
1847	PARAMETER (kstable=0.002)
1848	ccc PARAMETER (kstable=0.001)
1849	REAL mixlen ! constante controlant longueur de melange
1850	PARAMETER (mixlen=35.0)
1851	REAL seuil ! au-dela l'inversion est consideree trop faible
1852	PARAMETER (seuil=-0.02)
1853	ccc PARAMETER (seuil=-0.04)
1854	ccc PARAMETER (seuil=-0.06)
1855	ccc PARAMETER (seuil=-0.09)
1856	c
1857	c Variables locales:
1858	c
1859	INTEGER i, k, invb(knon)
1860	REAL zl2(knon)
1861	REAL zdthmin(knon), zdthdp
1862	c
1863	c Initialiser les sorties
1864	c
1865	DO k = 1, klev
1866	DO i = 1, knon
1867	pcfm(i,k) = 0.0
1868	pcfh(i,k) = 0.0
1869	ENDDO
1870	ENDDO
1871	c
1872	c Chercher la zone d'inversion forte
1873	c
1874	DO i = 1, knon
1875	invb(i) = klev
1876	zdthmin(i)=0.0
1877	ENDDO
1878	DO k = 2, klev/2-1
1879	DO i = 1, knon
1880	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1881	. - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
1882	zdthdp = zdthdp * 100.0
1883	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1884	. zdthdp.LT.zdthmin(i) ) THEN
1885	zdthmin(i) = zdthdp
1886	invb(i) = k
1887	ENDIF
1888	ENDDO
1889	ENDDO
1890	c
1891	c Introduire une diffusion:
1892	c
1893	DO k = 2, klev
1894	DO i = 1, knon
1895	cIM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean
1896	cIM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
1897	!IM cf JLD/ GKtest TERkz2
1898	! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre
1899	! fin GKtest
1900	IF ( (nsrf.EQ.is_oce) .AND. ! si on est sur ocean et si
1901	. ( (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
1902	. (zdthmin(i).GT.seuil) ) )THEN ! si l'inversion est trop faible
1903	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1))
1904	. /(paprs(i,2)-paprs(i,klev+1)) ))**2
1905	pcfm(i,k)= zl2(i)* kstable
1906	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1907	ENDIF
1908	ENDDO
1909	ENDDO
1910	c
1911	RETURN
1912	END
1913	SUBROUTINE calbeta(dtime,indice,knon,snow,qsol,
1914	. vbeta,vcal,vdif)
1915	IMPLICIT none
1916	c======================================================================
1917	c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
1918	c date: 19940414
1919	c======================================================================
1920	c
1921	c Calculer quelques parametres pour appliquer la couche limite
1922	c ------------------------------------------------------------
1923	#include "dimensions.h"
1924	#include "dimphy.h"
1925	#include "YOMCST.h"
1926	#include "indicesol.h"
1927	#include "iniprint.h"
1928	REAL tau_gl ! temps de relaxation pour la glace de mer
1929	ccc PARAMETER (tau_gl=86400.0*30.0)
1930	PARAMETER (tau_gl=86400.0*5.0)
1931	REAL mx_eau_sol
1932	PARAMETER (mx_eau_sol=150.0)
1933	c
1934	REAL calsol, calsno, calice ! epaisseur du sol: 0.15 m
1935	PARAMETER (calsol=1.0/(2.5578E+06*0.15))
1936	PARAMETER (calsno=1.0/(2.3867E+06*0.15))
1937	PARAMETER (calice=1.0/(5.1444E+06*0.15))
1938	C
1939	INTEGER i
1940	c
1941	REAL dtime
1942	REAL snow(klon), qsol(klon)
1943	INTEGER indice, knon
1944	C
1945	REAL vbeta(klon)
1946	REAL vcal(klon)
1947	REAL vdif(klon)
1948	C
1949
1950	IF (indice.EQ.is_oce) THEN
1951	DO i = 1, knon
1952	vcal(i) = 0.0
1953	vbeta(i) = 1.0
1954	vdif(i) = 0.0
1955	ENDDO
1956	ENDIF
1957	c
1958	IF (indice.EQ.is_sic) THEN
1959	DO i = 1, knon
1960	vcal(i) = calice
1961	IF (snow(i) .GT. 0.0) vcal(i) = calsno
1962	vbeta(i) = 1.0
1963	vdif(i) = 1.0/tau_gl
1964	ccc vdif(i) = calice/tau_gl ! c'etait une erreur
1965	ENDDO
1966	ENDIF
1967	c
1968	IF (indice.EQ.is_ter) THEN
1969	DO i = 1, knon
1970	vcal(i) = calsol
1971	IF (snow(i) .GT. 0.0) vcal(i) = calsno
1972	vbeta(i) = MIN(2.0*qsol(i)/mx_eau_sol, 1.0)
1973	vdif(i) = 0.0
1974	ENDDO
1975	ENDIF
1976	c
1977	IF (indice.EQ.is_lic) THEN
1978	DO i = 1, knon
1979	vcal(i) = calice
1980	IF (snow(i) .GT. 0.0) vcal(i) = calsno
1981	vbeta(i) = 1.0
1982	vdif(i) = 0.0
1983	ENDDO
1984	ENDIF
1985	c
1986	RETURN
1987	END
1988	C======================================================================
1989	SUBROUTINE nonlocal(knon, paprs, pplay,
1990	. tsol,beta,u,v,t,q,
1991	. cd_h, cd_m, pcfh, pcfm, cgh, cgq)
1992	IMPLICIT none
1993	c======================================================================
1994	c Laurent Li (LMD/CNRS), le 30 septembre 1998
1995	c Couche limite non-locale. Adaptation du code du CCM3.
1996	c Code non teste, donc a ne pas utiliser.
1997	c======================================================================
1998	c Nonlocal scheme that determines eddy diffusivities based on a
1999	c diagnosed boundary layer height and a turbulent velocity scale.
2000	c Also countergradient effects for heat and moisture are included.
2001	c
2002	c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
2003	c Local versus nonlocal boundary-layer diffusion in a global climate
2004	c model. J. of Climate, vol. 6, 1825-1842.
2005	c======================================================================
2006	#include "dimensions.h"
2007	#include "dimphy.h"
2008	#include "YOMCST.h"
2009	#include "iniprint.h"
2010	c
2011	c Arguments:
2012	c
2013	INTEGER knon ! nombre de points a calculer
2014	REAL tsol(klon) ! temperature du sol (K)
2015	REAL beta(klon) ! efficacite d'evaporation (entre 0 et 1)
2016	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
2017	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
2018	REAL u(klon,klev) ! vitesse U (m/s)
2019	REAL v(klon,klev) ! vitesse V (m/s)
2020	REAL t(klon,klev) ! temperature (K)
2021	REAL q(klon,klev) ! vapeur d'eau (kg/kg)
2022	REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
2023	REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
2024	c
2025	INTEGER isommet
2026	PARAMETER (isommet=klev)
2027	REAL vk
2028	PARAMETER (vk=0.40)
2029	REAL ricr
2030	PARAMETER (ricr=0.4)
2031	REAL fak
2032	PARAMETER (fak=8.5)
2033	REAL fakn
2034	PARAMETER (fakn=7.2)
2035	REAL onet
2036	PARAMETER (onet=1.0/3.0)
2037	REAL t_coup
2038	PARAMETER(t_coup=273.15)
2039	REAL zkmin
2040	PARAMETER (zkmin=0.01)
2041	REAL betam
2042	PARAMETER (betam=15.0)
2043	REAL betah
2044	PARAMETER (betah=15.0)
2045	REAL betas
2046	PARAMETER (betas=5.0)
2047	REAL sffrac
2048	PARAMETER (sffrac=0.1)
2049	REAL binm
2050	PARAMETER (binm=betam*sffrac)
2051	REAL binh
2052	PARAMETER (binh=betah*sffrac)
2053	REAL ccon
2054	PARAMETER (ccon=faksffracvk)
2055	c
2056	REAL z(klon,klev)
2057	REAL pcfm(klon,klev), pcfh(klon,klev)
2058	c
2059	INTEGER i, k
2060	REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
2061	REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
2062	REAL khfs(klon) ! surface kinematic heat flux [mK/s]
2063	REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
2064	REAL heatv(klon) ! surface virtual heat flux
2065	REAL ustar(klon)
2066	REAL rino(klon,klev) ! bulk Richardon no. from level to ref lev
2067	LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
2068	LOGICAL stblev(klon) ! stable pbl with levels within pbl
2069	LOGICAL unslev(klon) ! unstbl pbl with levels within pbl
2070	LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr
2071	LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr
2072	LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
2073	REAL pblh(klon)
2074	REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
2075	REAL cgq(klon,2:klev) ! counter-gradient term for constituents
2076	REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
2077	REAL obklen(klon)
2078	REAL ztvd, ztvu, zdu2
2079	REAL therm(klon) ! thermal virtual temperature excess
2080	REAL phiminv(klon) ! inverse phi function for momentum
2081	REAL phihinv(klon) ! inverse phi function for heat
2082	REAL wm(klon) ! turbulent velocity scale for momentum
2083	REAL fak1(klon) ! kustarpblh
2084	REAL fak2(klon) ! kwmpblh
2085	REAL fak3(klon) ! fakn*wstr/wm
2086	REAL pblk(klon) ! level eddy diffusivity for momentum
2087	REAL pr(klon) ! Prandtl number for eddy diffusivities
2088	REAL zl(klon) ! zmzp / Obukhov length
2089	REAL zh(klon) ! zmzp / pblh
2090	REAL zzh(klon) ! (1-(zmzp/pblh))**2
2091	REAL wstr(klon) ! w*, convective velocity scale
2092	REAL zm(klon) ! current level height
2093	REAL zp(klon) ! current level height + one level up
2094	REAL zcor, zdelta, zcvm5, zxqs
2095	REAL fac, pblmin, zmzp, term
2096	c
2097	#include "YOETHF.h"
2098	#include "FCTTRE.h"
2099	c
2100	c Initialisation
2101	c
2102	DO i = 1, klon
2103	pcfh(i,1) = cd_h(i)
2104	pcfm(i,1) = cd_m(i)
2105	ENDDO
2106	DO k = 2, klev
2107	DO i = 1, klon
2108	pcfh(i,k) = zkmin
2109	pcfm(i,k) = zkmin
2110	cgs(i,k) = 0.0
2111	cgh(i,k) = 0.0
2112	cgq(i,k) = 0.0
2113	ENDDO
2114	ENDDO
2115	c
2116	c Calculer les hauteurs de chaque couche
2117	c
2118	DO i = 1, knon
2119	z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
2120	. * (paprs(i,1)-pplay(i,1)) / RG
2121	ENDDO
2122	DO k = 2, klev
2123	DO i = 1, knon
2124	z(i,k) = z(i,k-1)
2125	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
2126	. * (pplay(i,k-1)-pplay(i,k)) / RG
2127	ENDDO
2128	ENDDO
2129	c
2130	DO i = 1, knon
2131	IF (thermcep) THEN
2132	zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
2133	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
2134	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
2135	zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
2136	zxqs=MIN(0.5,zxqs)
2137	zcor=1./(1.-retv*zxqs)
2138	zxqs=zxqs*zcor
2139	ELSE
2140	IF (tsol(i).LT.t_coup) THEN
2141	zxqs = qsats(tsol(i)) / paprs(i,1)
2142	ELSE
2143	zxqs = qsatl(tsol(i)) / paprs(i,1)
2144	ENDIF
2145	ENDIF
2146	zx_alf1 = 1.0
2147	zx_alf2 = 1.0 - zx_alf1
2148	zxt = (t(i,1)+z(i,1)RG/RCPD/(1.+RVTMP2q(i,1)))
2149	. (1.+RETVq(i,1))*zx_alf1
2150	. + (t(i,2)+z(i,2)RG/RCPD/(1.+RVTMP2q(i,2)))
2151	. (1.+RETVq(i,2))*zx_alf2
2152	zxu = u(i,1)zx_alf1+u(i,2)zx_alf2
2153	zxv = v(i,1)zx_alf1+v(i,2)zx_alf2
2154	zxq = q(i,1)zx_alf1+q(i,2)zx_alf2
2155	zxmod = 1.0+SQRT(zxu2+zxv2)
2156	khfs(i) = (tsol(i)(1.+RETVq(i,1))-zxt) zxmodcd_h(i)
2157	kqfs(i) = (zxqs-zxq) zxmodcd_h(i) * beta(i)
2158	heatv(i) = khfs(i) + 0.61zxtkqfs(i)
2159	taux = zxu zxmodcd_m(i)
2160	tauy = zxv zxmodcd_m(i)
2161	ustar(i) = SQRT(taux2+tauy2)
2162	ustar(i) = MAX(SQRT(ustar(i)),0.01)
2163	ENDDO
2164	c
2165	DO i = 1, knon
2166	rino(i,1) = 0.0
2167	check(i) = .TRUE.
2168	pblh(i) = z(i,1)
2169	obklen(i) = -t(i,1)ustar(i)3/(RGvk*heatv(i))
2170	ENDDO
2171
2172	C
2173	C PBL height calculation:
2174	C Search for level of pbl. Scan upward until the Richardson number between
2175	C the first level and the current level exceeds the "critical" value.
2176	C
2177	fac = 100.0
2178	DO k = 1, isommet
2179	DO i = 1, knon
2180	IF (check(i)) THEN
2181	zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
2182	zdu2 = max(zdu2,1.0e-20)
2183	ztvd =(t(i,k)+z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,k)))
2184	. (1.+RETVq(i,k))
2185	ztvu =(t(i,1)-z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
2186	. (1.+RETVq(i,1))
2187	rino(i,k) = (z(i,k)-z(i,1))RG(ztvd-ztvu)
2188	. /(zdu20.5(ztvd+ztvu))
2189	IF (rino(i,k).GE.ricr) THEN
2190	pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
2191	. (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
2192	check(i) = .FALSE.
2193	ENDIF
2194	ENDIF
2195	ENDDO
2196	ENDDO
2197
2198	C
2199	C Set pbl height to maximum value where computation exceeds number of
2200	C layers allowed
2201	C
2202	DO i = 1, knon
2203	if (check(i)) pblh(i) = z(i,isommet)
2204	ENDDO
2205	C
2206	C Improve estimate of pbl height for the unstable points.
2207	C Find unstable points (sensible heat flux is upward):
2208	C
2209	DO i = 1, knon
2210	IF (heatv(i) .GT. 0.) THEN
2211	unstbl(i) = .TRUE.
2212	check(i) = .TRUE.
2213	ELSE
2214	unstbl(i) = .FALSE.
2215	check(i) = .FALSE.
2216	ENDIF
2217	ENDDO
2218	C
2219	C For the unstable case, compute velocity scale and the
2220	C convective temperature excess:
2221	C
2222	DO i = 1, knon
2223	IF (check(i)) THEN
2224	phiminv(i) = (1.-binmpblh(i)/obklen(i))*onet
2225	wm(i)= ustar(i)*phiminv(i)
2226	therm(i) = heatv(i)*fak/wm(i)
2227	rino(i,1) = 0.0
2228	ENDIF
2229	ENDDO
2230	C
2231	C Improve pblh estimate for unstable conditions using the
2232	C convective temperature excess:
2233	C
2234	DO k = 1, isommet
2235	DO i = 1, knon
2236	IF (check(i)) THEN
2237	zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
2238	zdu2 = max(zdu2,1.0e-20)
2239	ztvd =(t(i,k)+z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,k)))
2240	. (1.+RETVq(i,k))
2241	ztvu =(t(i,1)+therm(i)-z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
2242	. (1.+RETVq(i,1))
2243	rino(i,k) = (z(i,k)-z(i,1))RG(ztvd-ztvu)
2244	. /(zdu20.5(ztvd+ztvu))
2245	IF (rino(i,k).GE.ricr) THEN
2246	pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
2247	. (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
2248	check(i) = .FALSE.
2249	ENDIF
2250	ENDIF
2251	ENDDO
2252	ENDDO
2253	C
2254	C Set pbl height to maximum value where computation exceeds number of
2255	C layers allowed
2256	C
2257	DO i = 1, knon
2258	if (check(i)) pblh(i) = z(i,isommet)
2259	ENDDO
2260	C
2261	C Points for which pblh exceeds number of pbl layers allowed;
2262	C set to maximum
2263	C
2264	DO i = 1, knon
2265	IF (check(i)) pblh(i) = z(i,isommet)
2266	ENDDO
2267	C
2268	C PBL height must be greater than some minimum mechanical mixing depth
2269	C Several investigators have proposed minimum mechanical mixing depth
2270	C relationships as a function of the local friction velocity, u*. We
2271	C make use of a linear relationship of the form h = c u* where c=700.
2272	C The scaling arguments that give rise to this relationship most often
2273	C represent the coefficient c as some constant over the local coriolis
2274	C parameter. Here we make use of the experimental results of Koracin
2275	C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
2276	C where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
2277	C latitude value for f so that c = 0.07/f = 700.
2278	C
2279	DO i = 1, knon
2280	pblmin = 700.0*ustar(i)
2281	pblh(i) = MAX(pblh(i),pblmin)
2282	ENDDO
2283	C
2284	C pblh is now available; do preparation for diffusivity calculation:
2285	C
2286	DO i = 1, knon
2287	pblk(i) = 0.0
2288	fak1(i) = ustar(i)pblh(i)vk
2289	C
2290	C Do additional preparation for unstable cases only, set temperature
2291	C and moisture perturbations depending on stability.
2292	C
2293	IF (unstbl(i)) THEN
2294	zxt=(t(i,1)-z(i,1)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
2295	. (1.+RETVq(i,1))
2296	phiminv(i) = (1. - binmpblh(i)/obklen(i))*onet
2297	phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
2298	wm(i) = ustar(i)*phiminv(i)
2299	fak2(i) = wm(i)pblh(i)vk
2300	wstr(i) = (heatv(i)RGpblh(i)/zxt)**onet
2301	fak3(i) = fakn*wstr(i)/wm(i)
2302	ENDIF
2303	ENDDO
2304
2305	C Main level loop to compute the diffusivities and
2306	C counter-gradient terms:
2307	C
2308	DO 1000 k = 2, isommet
2309	C
2310	C Find levels within boundary layer:
2311	C
2312	DO i = 1, knon
2313	unslev(i) = .FALSE.
2314	stblev(i) = .FALSE.
2315	zm(i) = z(i,k-1)
2316	zp(i) = z(i,k)
2317	IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
2318	IF (zm(i) .LT. pblh(i)) THEN
2319	zmzp = 0.5*(zm(i) + zp(i))
2320	zh(i) = zmzp/pblh(i)
2321	zl(i) = zmzp/obklen(i)
2322	zzh(i) = 0.
2323	IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
2324	C
2325	C stblev for points zm < plbh and stable and neutral
2326	C unslev for points zm < plbh and unstable
2327	C
2328	IF (unstbl(i)) THEN
2329	unslev(i) = .TRUE.
2330	ELSE
2331	stblev(i) = .TRUE.
2332	ENDIF
2333	ENDIF
2334	ENDDO
2335	C
2336	C Stable and neutral points; set diffusivities; counter-gradient
2337	C terms zero for stable case:
2338	C
2339	DO i = 1, knon
2340	IF (stblev(i)) THEN
2341	IF (zl(i).LE.1.) THEN
2342	pblk(i) = fak1(i)zh(i)zzh(i)/(1. + betas*zl(i))
2343	ELSE
2344	pblk(i) = fak1(i)zh(i)zzh(i)/(betas + zl(i))
2345	ENDIF
2346	pcfm(i,k) = pblk(i)
2347	pcfh(i,k) = pcfm(i,k)
2348	ENDIF
2349	ENDDO
2350	C
2351	C unssrf, unstable within surface layer of pbl
2352	C unsout, unstable within outer layer of pbl
2353	C
2354	DO i = 1, knon
2355	unssrf(i) = .FALSE.
2356	unsout(i) = .FALSE.
2357	IF (unslev(i)) THEN
2358	IF (zh(i).lt.sffrac) THEN
2359	unssrf(i) = .TRUE.
2360	ELSE
2361	unsout(i) = .TRUE.
2362	ENDIF
2363	ENDIF
2364	ENDDO
2365	C
2366	C Unstable for surface layer; counter-gradient terms zero
2367	C
2368	DO i = 1, knon
2369	IF (unssrf(i)) THEN
2370	term = (1. - betamzl(i))*onet
2371	pblk(i) = fak1(i)zh(i)zzh(i)*term
2372	pr(i) = term/sqrt(1. - betah*zl(i))
2373	ENDIF
2374	ENDDO
2375	C
2376	C Unstable for outer layer; counter-gradient terms non-zero:
2377	C
2378	DO i = 1, knon
2379	IF (unsout(i)) THEN
2380	pblk(i) = fak2(i)zh(i)zzh(i)
2381	cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
2382	cgh(i,k) = khfs(i)*cgs(i,k)
2383	pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
2384	cgq(i,k) = kqfs(i)*cgs(i,k)
2385	ENDIF
2386	ENDDO
2387	C
2388	C For all unstable layers, set diffusivities
2389	C
2390	DO i = 1, knon
2391	IF (unslev(i)) THEN
2392	pcfm(i,k) = pblk(i)
2393	pcfh(i,k) = pblk(i)/pr(i)
2394	ENDIF
2395	ENDDO
2396	1000 continue ! end of level loop
2397
2398	RETURN
2399	END

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

Download in other formats:

Original Format