Context Navigation

clmain.F @ 3733

Last change on this file since 3733 was 644, checked in by Laurent Fairhead, 19 years ago
Synchronisation avec tous les diagnostiques de Ionela IM Inclusion du slab ocean IM 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	fqcalving(i,nsrf) = y_fqcalving(j)
791	ffonte(i,nsrf) = y_ffonte(j)
792	cdragh(i) = cdragh(i) + ycoefh(j,1)
793	cdragm(i) = cdragm(i) + ycoefm(j,1)
794	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
795	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
796	zu1(i) = zu1(i) + yu1(j)
797	zv1(i) = zv1(i) + yv1(j)
798	END DO
799	IF ( nsrf .eq. is_ter ) THEN
800	DO j = 1, knon
801	i = ni(j)
802	qsol(i) = yqsol(j)
803	END DO
804	END IF
805	IF ( nsrf .eq. is_lic ) THEN
806	DO j = 1, knon
807	i = ni(j)
808	run_off_lic_0(i) = y_run_off_lic_0(j)
809	END DO
810	END IF
811	c$$$ PB ajout pour soil
812	ftsoil(:,:,nsrf) = 0.
813	DO k = 1, nsoilmx
814	DO j = 1, knon
815	i = ni(j)
816	ftsoil(i, k, nsrf) = ytsoil(j,k)
817	END DO
818	END DO
819	c
820	#ifdef CRAY
821	DO k = 1, klev
822	DO j = 1, knon
823	i = ni(j)
824	#else
825	DO j = 1, knon
826	i = ni(j)
827	DO k = 1, klev
828	#endif
829	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
830	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
831	c$$$ PB flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
832	c$$$ flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
833	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
834	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
835	c$$$ PB flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
836	c$$$ flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
837	zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
838	ENDDO
839	ENDDO
840	c
841	c
842	#undef T2m
843	#define T2m
844	#ifdef T2m
845	ccc diagnostic t,q a 2m et u, v a 10m
846	c
847	DO j=1, knon
848	i = ni(j)
849	uzon(j) = yu(j,1) + y_d_u(j,1)
850	vmer(j) = yv(j,1) + y_d_v(j,1)
851	tair1(j) = yt(j,1) + y_d_t(j,1)
852	qair1(j) = yq(j,1) + y_d_q(j,1)
853	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1)))
854	& * (ypaprs(j,1)-ypplay(j,1))
855	tairsol(j) = yts(j) + y_d_ts(j)
856	rugo1(j) = yrugos(j)
857	IF(nsrf.EQ.is_oce) THEN
858	rugo1(j) = rugos(i,nsrf)
859	ENDIF
860	psfce(j)=ypaprs(j,1)
861	patm(j)=ypplay(j,1)
862	c
863	qairsol(j) = yqsurf(j)
864	ENDDO
865	c
866	CALL stdlevvar(klon, knon, nsrf, zxli,
867	& uzon, vmer, tair1, qair1, zgeo1,
868	& tairsol, qairsol, rugo1, psfce, patm,
869	cIM & yt2m, yq2m, yu10m)
870	& yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
871	cIM 081204 END
872	c
873	c
874	DO j=1, knon
875	i = ni(j)
876	t2m(i,nsrf)=yt2m(j)
877
878	c
879	q2m(i,nsrf)=yq2m(j)
880	c
881	c u10m, v10m : composantes du vent a 10m sans spirale de Ekman
882	u10m(i,nsrf)=(yu10m(j) * uzon(j))/sqrt(uzon(j)2+vmer(j)2)
883	v10m(i,nsrf)=(yu10m(j) * vmer(j))/sqrt(uzon(j)2+vmer(j)2)
884	c
885	ENDDO
886	c
887	cIM cf AM : pbl, HBTM
888	DO i = 1, knon
889	y_cd_h(i) = ycoefh(i,1)
890	y_cd_m(i) = ycoefm(i,1)
891	ENDDO
892	c print*,'appel hbtm2'
893	CALL HBTM(knon, ypaprs, ypplay,
894	. yt2m,yt10m,yq2m,yq10m,yustar,
895	. y_flux_t,y_flux_q,yu,yv,yt,yq,
896	. ypblh,ycapCL,yoliqCL,ycteiCL,ypblT,
897	. ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
898	c print*,'fin hbtm2'
899	c
900	DO j=1, knon
901	i = ni(j)
902	pblh(i,nsrf) = ypblh(j)
903	plcl(i,nsrf) = ylcl(j)
904	capCL(i,nsrf) = ycapCL(j)
905	oliqCL(i,nsrf) = yoliqCL(j)
906	cteiCL(i,nsrf) = ycteiCL(j)
907	pblT(i,nsrf) = ypblT(j)
908	therm(i,nsrf) = ytherm(j)
909	trmb1(i,nsrf) = ytrmb1(j)
910	trmb2(i,nsrf) = ytrmb2(j)
911	trmb3(i,nsrf) = ytrmb3(j)
912	ENDDO
913	c
914	#else
915	cIM cf AM : pbl, HBTM
916	DO j=1, knon
917	i = ni(j)
918	pblh(i,nsrf) = 0.
919	plcl(i,nsrf) = 0.
920	capCL(i,nsrf) = 0.
921	oliqCL(i,nsrf) = 0.
922	cteiCL(i,nsrf) = 0.
923	pblT(i,nsrf) = 0.
924	therm(i,nsrf) = 0.
925	trmb1(i,nsrf) = 0.
926	trmb2(i,nsrf) = 0.
927	trmb3(i,nsrf) = 0.
928	ENDDO
929	DO j=1, knon
930	i = ni(j)
931	t2m(i,nsrf)=0.
932	q2m(i,nsrf)=0.
933	u10m(i,nsrf)=0.
934	v10m(i,nsrf)=0.
935	ENDDO
936	#endif
937
938	do j=1,knon
939	do k=1,klev+1
940	i=ni(j)
941	q2(i,k,nsrf)=yq2(j,k)
942	enddo
943	enddo
944	cIM "slab" ocean
945	IF(OCEAN.EQ.'slab '.OR.OCEAN.EQ.'force ') THEN
946	IF (nsrf.EQ.is_oce) THEN
947	DO j = 1, knon
948	c on projette sur la grille globale
949	i = ni(j)
950	IF(pctsrf_new(i,is_oce).GT.epsfra) THEN
951	flux_o(i) = y_flux_o(j)
952	ELSE
953	flux_o(i) = 0.
954	ENDIF
955	ENDDO
956	ENDIF
957	c
958	IF (nsrf.EQ.is_sic) THEN
959	DO j = 1, knon
960	i = ni(j)
961	cIM 230604 on pondere lorsque l'on fait le bilan au sol : flux_g(i) = y_flux_g(j)*ypct(j)
962	IF(pctsrf_new(i,is_sic).GT.epsfra) THEN
963	flux_g(i) = y_flux_g(j)
964	ELSE
965	flux_g(i) = 0.
966	ENDIF
967	ENDDO
968	ENDIF !nsrf.EQ.is_sic
969	ENDIF !OCEAN
970	c
971	IF(OCEAN.EQ.'slab ') THEN
972	IF(nsrf.EQ.is_oce) then
973	tslab(1:klon) = ytslab(1:klon)
974	seaice(1:klon) = y_seaice(1:klon)
975	ENDIF !nsrf
976	ENDIF !OCEAN
977	99999 CONTINUE
978	C
979	C On utilise les nouvelles surfaces
980	C A rajouter: conservation de l'albedo
981	C
982	rugos(:,is_oce) = rugmer
983	pctsrf = pctsrf_new
984
985	RETURN
986	END
987	SUBROUTINE clqh(dtime,itime, date0,jour,debut,lafin,
988	e rlon, rlat, cufi, cvfi,
989	e knon, nisurf, knindex, pctsrf,
990	$ soil_model,tsoil,qsol,
991	e ok_veget, ocean, npas, nexca,
992	e rmu0, co2_ppm, rugos, rugoro,
993	e u1lay,v1lay,coef,
994	e t,q,ts,paprs,pplay,
995	e delp,radsol,albedo,alblw,snow,qsurf,
996	e precip_rain, precip_snow, fder, taux, tauy,
997	c -- LOOP
998	e ywindsp,
999	c -- LOOP
1000	$ sollw, sollwdown, swnet,fluxlat,
1001	s pctsrf_new, agesno,
1002	s d_t, d_q, d_ts, z0_new,
1003	s flux_t, flux_q,dflux_s,dflux_l,
1004	s fqcalving,ffonte,run_off_lic_0,
1005	cIM "slab" ocean
1006	s flux_o,flux_g,tslab,seaice)
1007
1008	USE interface_surf
1009
1010	IMPLICIT none
1011	c======================================================================
1012	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
1013	c Objet: diffusion verticale de "q" et de "h"
1014	c======================================================================
1015	#include "dimensions.h"
1016	#include "dimphy.h"
1017	#include "YOMCST.h"
1018	#include "YOETHF.h"
1019	#include "FCTTRE.h"
1020	#include "indicesol.h"
1021	#include "dimsoil.h"
1022	#include "iniprint.h"
1023
1024	c Arguments:
1025	INTEGER knon
1026	REAL dtime ! intervalle du temps (s)
1027	real date0
1028	REAL u1lay(klon) ! vitesse u de la 1ere couche (m/s)
1029	REAL v1lay(klon) ! vitesse v de la 1ere couche (m/s)
1030	REAL coef(klon,klev) ! le coefficient d'echange (m**2/s)
1031	c multiplie par le cisaillement du
1032	c vent (dV/dz); la premiere valeur
1033	c indique la valeur de Cdrag (sans unite)
1034	REAL t(klon,klev) ! temperature (K)
1035	REAL q(klon,klev) ! humidite specifique (kg/kg)
1036	REAL ts(klon) ! temperature du sol (K)
1037	REAL evap(klon) ! evaporation au sol
1038	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
1039	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
1040	REAL delp(klon,klev) ! epaisseur de couche en pression (Pa)
1041	REAL radsol(klon) ! ray. net au sol (Solaire+IR) W/m2
1042	REAL albedo(klon) ! albedo de la surface
1043	REAL alblw(klon)
1044	REAL snow(klon) ! hauteur de neige
1045	REAL qsurf(klon) ! humidite de l'air au dessus de la surface
1046	real precip_rain(klon), precip_snow(klon)
1047	REAL agesno(klon)
1048	REAL rugoro(klon)
1049	REAL run_off_lic_0(klon)! runof glacier au pas de temps precedent
1050	integer jour ! jour de l'annee en cours
1051	real rmu0(klon) ! cosinus de l'angle solaire zenithal
1052	real rugos(klon) ! rugosite
1053	integer knindex(klon)
1054	real pctsrf(klon,nbsrf)
1055	real rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
1056	logical ok_veget
1057	REAL co2_ppm ! taux CO2 atmosphere
1058	character*6 ocean
1059	integer npas, nexca
1060	c -- LOOP
1061	REAL yu10mx(klon)
1062	REAL yu10my(klon)
1063	REAL ywindsp(klon)
1064	c -- LOOP
1065
1066
1067	c
1068	REAL d_t(klon,klev) ! incrementation de "t"
1069	REAL d_q(klon,klev) ! incrementation de "q"
1070	REAL d_ts(klon) ! incrementation de "ts"
1071	REAL flux_t(klon,klev) ! (diagnostic) flux de la chaleur
1072	c sensible, flux de Cp*T, positif vers
1073	c le bas: j/(m**2 s) c.a.d.: W/m2
1074	REAL flux_q(klon,klev) ! flux de la vapeur d'eau:kg/(m**2 s)
1075	REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
1076	REAL dflux_l(klon) ! derivee du flux latent dF/dTs
1077	cIM cf JLD
1078	c Flux thermique utiliser pour fondre la neige
1079	REAL ffonte(klon)
1080	c Flux d'eau "perdue" par la surface et nécessaire pour que limiter la
1081	c hauteur de neige, en kg/m2/s
1082	REAL fqcalving(klon)
1083	cIM "slab" ocean
1084	REAL tslab(klon) !temperature du slab ocean (K) (OCEAN='slab ')
1085	REAL seaice(klon) ! glace de mer en kg/m2
1086	REAL flux_o(klon) ! flux entre l'ocean et l'atmosphere W/m2
1087	REAL flux_g(klon) ! flux entre l'ocean et la glace de mer W/m2
1088	c
1089	c======================================================================
1090	REAL t_grnd ! temperature de rappel pour glace de mer
1091	PARAMETER (t_grnd=271.35)
1092	REAL t_coup
1093	PARAMETER(t_coup=273.15)
1094	c======================================================================
1095	INTEGER i, k
1096	REAL zx_cq(klon,klev)
1097	REAL zx_dq(klon,klev)
1098	REAL zx_ch(klon,klev)
1099	REAL zx_dh(klon,klev)
1100	REAL zx_buf1(klon)
1101	REAL zx_buf2(klon)
1102	REAL zx_coef(klon,klev)
1103	REAL local_h(klon,klev) ! enthalpie potentielle
1104	REAL local_q(klon,klev)
1105	REAL local_ts(klon)
1106	REAL psref(klon) ! pression de reference pour temperature potent.
1107	REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
1108	c======================================================================
1109	c contre-gradient pour la vapeur d'eau: (kg/kg)/metre
1110	REAL gamq(klon,2:klev)
1111	c contre-gradient pour la chaleur sensible: Kelvin/metre
1112	REAL gamt(klon,2:klev)
1113	REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)
1114	REAL zdelz
1115	c======================================================================
1116	#include "compbl.h"
1117	c======================================================================
1118	c Rajout pour l'interface
1119	integer itime
1120	integer nisurf
1121	logical debut, lafin
1122	real zlev1(klon)
1123	real fder(klon), taux(klon), tauy(klon)
1124	real temp_air(klon), spechum(klon)
1125	real epot_air(klon), ccanopy(klon)
1126	real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
1127	real petBcoef(klon), peqBcoef(klon)
1128	real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
1129	real p1lay(klon)
1130	c$$$C PB ajout pour soil
1131	LOGICAL soil_model
1132	REAL tsoil(klon, nsoilmx)
1133	REAL qsol(klon)
1134
1135	! Parametres de sortie
1136	real fluxsens(klon), fluxlat(klon)
1137	real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
1138	real emis_new(klon), z0_new(klon)
1139	real pctsrf_new(klon,nbsrf)
1140	c JLD
1141	real zzpk
1142	C
1143	character (len = 20) :: modname = 'Debut clqh'
1144	LOGICAL check
1145	PARAMETER (check=.false.)
1146	C
1147	if (check) THEN
1148	write(,) modname,' nisurf=',nisurf
1149	call flush(6)
1150	endif
1151	c
1152	if (check) THEN
1153	WRITE(,)' qsurf (min, max)'
1154	$ , minval(qsurf(1:knon)), maxval(qsurf(1:knon))
1155	call flush(6)
1156	ENDIF
1157	C
1158	C
1159	if (iflag_pbl.eq.1) then
1160	do k = 3, klev
1161	do i = 1, knon
1162	gamq(i,k)= 0.0
1163	gamt(i,k)= -1.0e-03
1164	enddo
1165	enddo
1166	do i = 1, knon
1167	gamq(i,2) = 0.0
1168	gamt(i,2) = -2.5e-03
1169	enddo
1170	else
1171	do k = 2, klev
1172	do i = 1, knon
1173	gamq(i,k) = 0.0
1174	gamt(i,k) = 0.0
1175	enddo
1176	enddo
1177	endif
1178
1179	DO i = 1, knon
1180	psref(i) = paprs(i,1) !pression de reference est celle au sol
1181	local_ts(i) = ts(i)
1182	ENDDO
1183	DO k = 1, klev
1184	DO i = 1, knon
1185	zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
1186	zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
1187	local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)
1188	local_q(i,k) = q(i,k)
1189	ENDDO
1190	ENDDO
1191	c
1192	c Convertir les coefficients en variables convenables au calcul:
1193	c
1194	c
1195	DO k = 2, klev
1196	DO i = 1, knon
1197	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
1198	. (paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
1199	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
1200	ENDDO
1201	ENDDO
1202	c
1203	c Preparer les flux lies aux contre-gardients
1204	c
1205	DO k = 2, klev
1206	DO i = 1, knon
1207	zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k)
1208	. *(pplay(i,k-1)-pplay(i,k))
1209	z_gamaq(i,k) = gamq(i,k) * zdelz
1210	z_gamah(i,k) = gamt(i,k) * zdelz RCPD zx_pkh(i,k)
1211	ENDDO
1212	ENDDO
1213	DO i = 1, knon
1214	zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
1215	zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev)
1216	. -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)
1217	zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
1218	c
1219	zzpk=(pplay(i,klev)/psref(i))**RKAPPA
1220	zx_buf2(i) = zzpk*delp(i,klev) + zx_coef(i,klev)
1221	zx_ch(i,klev) = (local_h(i,klev)zzpkdelp(i,klev)
1222	. -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
1223	zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
1224	ENDDO
1225	DO k = klev-1, 2 , -1
1226	DO i = 1, knon
1227	zx_buf1(i) = delp(i,k)+zx_coef(i,k)
1228	. +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))
1229	zx_cq(i,k) = (local_q(i,k)*delp(i,k)
1230	. +zx_coef(i,k+1)*zx_cq(i,k+1)
1231	. +zx_coef(i,k+1)*z_gamaq(i,k+1)
1232	. -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)
1233	zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
1234	c
1235	zzpk=(pplay(i,k)/psref(i))**RKAPPA
1236	zx_buf2(i) = zzpk*delp(i,k)+zx_coef(i,k)
1237	. +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
1238	zx_ch(i,k) = (local_h(i,k)zzpkdelp(i,k)
1239	. +zx_coef(i,k+1)*zx_ch(i,k+1)
1240	. +zx_coef(i,k+1)*z_gamah(i,k+1)
1241	. -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
1242	zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
1243	ENDDO
1244	ENDDO
1245	C
1246	C nouvelle formulation JL Dufresne
1247	C
1248	C q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt
1249	C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
1250	C
1251	DO i = 1, knon
1252	zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))
1253	zx_cq(i,1) = (local_q(i,1)*delp(i,1)
1254	. +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2)))
1255	. /zx_buf1(i)
1256	zx_dq(i,1) = -1. * RG / zx_buf1(i)
1257	c
1258	zzpk=(pplay(i,1)/psref(i))**RKAPPA
1259	zx_buf2(i) = zzpkdelp(i,1) + zx_coef(i,2)(1.-zx_dh(i,2))
1260	zx_ch(i,1) = (local_h(i,1)zzpkdelp(i,1)
1261	. +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
1262	. /zx_buf2(i)
1263	zx_dh(i,1) = -1. * RG / zx_buf2(i)
1264	ENDDO
1265
1266	C Appel a interfsurf (appel generique) routine d'interface avec la surface
1267
1268	c initialisation
1269	petAcoef =0.
1270	peqAcoef = 0.
1271	petBcoef =0.
1272	peqBcoef = 0.
1273	p1lay =0.
1274
1275	c do i = 1, knon
1276	petAcoef(1:knon) = zx_ch(1:knon,1)
1277	peqAcoef(1:knon) = zx_cq(1:knon,1)
1278	petBcoef(1:knon) = zx_dh(1:knon,1)
1279	peqBcoef(1:knon) = zx_dq(1:knon,1)
1280	tq_cdrag(1:knon) =coef(1:knon,1)
1281	temp_air(1:knon) =t(1:knon,1)
1282	epot_air(1:knon) =local_h(1:knon,1)
1283	spechum(1:knon)=q(1:knon,1)
1284	p1lay(1:knon) = pplay(1:knon,1)
1285	zlev1(1:knon) = delp(1:knon,1)
1286	c swnet = swdown * (1. - albedo)
1287	c
1288	cIM swdown=flux SW incident sur terres
1289	cIM swdown=flux SW net sur les autres surfaces
1290	cIM swdown(1:knon) = swnet(1:knon)
1291	if(nisurf.eq.is_ter) THEN
1292	swdown(1:knon) = swnet(1:knon)/(1-albedo(1:knon))
1293	else
1294	swdown(1:knon) = swnet(1:knon)
1295	endif
1296	c enddo
1297	ccanopy = co2_ppm
1298
1299	CALL interfsurf(itime, dtime, date0, jour, rmu0,
1300	e klon, iim, jjm, nisurf, knon, knindex, pctsrf,
1301	e rlon, rlat, cufi, cvfi,
1302	e debut, lafin, ok_veget, soil_model, nsoilmx,tsoil, qsol,
1303	e zlev1, u1lay, v1lay, temp_air, spechum, epot_air, ccanopy,
1304	e tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
1305	e precip_rain, precip_snow, sollw, sollwdown, swnet, swdown,
1306	e fder, taux, tauy,
1307	c -- LOOP
1308	e ywindsp,
1309	c -- LOOP
1310	e rugos, rugoro,
1311	e albedo, snow, qsurf,
1312	e ts, p1lay, psref, radsol,
1313	e ocean, npas, nexca, zmasq,
1314	s evap, fluxsens, fluxlat, dflux_l, dflux_s,
1315	s tsol_rad, tsurf_new, alb_new, alblw, emis_new, z0_new,
1316	s pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0,
1317	cIM "slab" ocean
1318	s flux_o, flux_g, tslab, seaice)
1319
1320
1321	do i = 1, knon
1322	flux_t(i,1) = fluxsens(i)
1323	flux_q(i,1) = - evap(i)
1324	d_ts(i) = tsurf_new(i) - ts(i)
1325	albedo(i) = alb_new(i)
1326	enddo
1327
1328	c==== une fois on a zx_h_ts, on peut faire l'iteration ========
1329	DO i = 1, knon
1330	local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)flux_t(i,1)dtime
1331	local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)flux_q(i,1)dtime
1332	ENDDO
1333	DO k = 2, klev
1334	DO i = 1, knon
1335	local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)
1336	local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
1337	ENDDO
1338	ENDDO
1339	c======================================================================
1340	c== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
1341	c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
1342	DO k = 2, klev
1343	DO i = 1, knon
1344	flux_q(i,k) = (zx_coef(i,k)/RG/dtime)
1345	. * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))
1346	flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
1347	. * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
1348	. / zx_pkh(i,k)
1349	ENDDO
1350	ENDDO
1351	c======================================================================
1352	C Calcul tendances
1353	DO k = 1, klev
1354	DO i = 1, knon
1355	d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)
1356	d_q(i,k) = local_q(i,k) - q(i,k)
1357	ENDDO
1358	ENDDO
1359	c
1360
1361	RETURN
1362	END
1363	SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
1364	e paprs,pplay,delp,
1365	s d_ven,flux_v)
1366	IMPLICIT none
1367	c======================================================================
1368	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
1369	c Objet: diffusion vertical de la vitesse "ven"
1370	c======================================================================
1371	c Arguments:
1372	c dtime----input-R- intervalle du temps (en second)
1373	c u1lay----input-R- vent u de la premiere couche (m/s)
1374	c v1lay----input-R- vent v de la premiere couche (m/s)
1375	c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
1376	c le cisaillement du vent (dV/dz); la premiere
1377	c valeur indique la valeur de Cdrag (sans unite)
1378	c t--------input-R- temperature (K)
1379	c ven------input-R- vitesse horizontale (m/s)
1380	c paprs----input-R- pression a inter-couche (Pa)
1381	c pplay----input-R- pression au milieu de couche (Pa)
1382	c delp-----input-R- epaisseur de couche (Pa)
1383	c
1384	c
1385	c d_ven----output-R- le changement de "ven"
1386	c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
1387	c======================================================================
1388	#include "dimensions.h"
1389	#include "dimphy.h"
1390	#include "iniprint.h"
1391	INTEGER knon
1392	REAL dtime
1393	REAL u1lay(klon), v1lay(klon)
1394	REAL coef(klon,klev)
1395	REAL t(klon,klev), ven(klon,klev)
1396	REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
1397	REAL d_ven(klon,klev)
1398	REAL flux_v(klon,klev)
1399	c======================================================================
1400	#include "YOMCST.h"
1401	c======================================================================
1402	INTEGER i, k
1403	REAL zx_cv(klon,2:klev)
1404	REAL zx_dv(klon,2:klev)
1405	REAL zx_buf(klon)
1406	REAL zx_coef(klon,klev)
1407	REAL local_ven(klon,klev)
1408	REAL zx_alf1(klon), zx_alf2(klon)
1409	c======================================================================
1410	DO k = 1, klev
1411	DO i = 1, knon
1412	local_ven(i,k) = ven(i,k)
1413	ENDDO
1414	ENDDO
1415	c======================================================================
1416	DO i = 1, knon
1417	ccc zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
1418	zx_alf1(i) = 1.0
1419	zx_alf2(i) = 1.0 - zx_alf1(i)
1420	zx_coef(i,1) = coef(i,1)
1421	. * (1.0+SQRT(u1lay(i)2+v1lay(i)2))
1422	. * pplay(i,1)/(RD*t(i,1))
1423	zx_coef(i,1) = zx_coef(i,1) * dtime*RG
1424	ENDDO
1425	c======================================================================
1426	DO k = 2, klev
1427	DO i = 1, knon
1428	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
1429	. (paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
1430	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
1431	ENDDO
1432	ENDDO
1433	c======================================================================
1434	DO i = 1, knon
1435	zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
1436	zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
1437	zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
1438	. /zx_buf(i)
1439	ENDDO
1440	DO k = 3, klev
1441	DO i = 1, knon
1442	zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
1443	. + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
1444	zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
1445	. +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
1446	zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
1447	ENDDO
1448	ENDDO
1449	DO i = 1, knon
1450	local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
1451	. +zx_coef(i,klev)*zx_cv(i,klev) )
1452	. / ( delp(i,klev) + zx_coef(i,klev)
1453	. -zx_coef(i,klev)*zx_dv(i,klev) )
1454	ENDDO
1455	DO k = klev-1, 1, -1
1456	DO i = 1, knon
1457	local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
1458	ENDDO
1459	ENDDO
1460	c======================================================================
1461	c== flux_v est le flux de moment angulaire (positif vers bas)
1462	c== dont l'unite est: (kg m/s)/(m**2 s)
1463	DO i = 1, knon
1464	flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
1465	. (local_ven(i,1)zx_alf1(i)
1466	. +local_ven(i,2)*zx_alf2(i))
1467	ENDDO
1468	DO k = 2, klev
1469	DO i = 1, knon
1470	flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
1471	. * (local_ven(i,k)-local_ven(i,k-1))
1472	ENDDO
1473	ENDDO
1474	c
1475	DO k = 1, klev
1476	DO i = 1, knon
1477	d_ven(i,k) = local_ven(i,k) - ven(i,k)
1478	ENDDO
1479	ENDDO
1480	c
1481	RETURN
1482	END
1483	SUBROUTINE coefkz(nsrf, knon, paprs, pplay,
1484	cIM 261103
1485	. ksta, ksta_ter,
1486	cIM 261103
1487	. ts, rugos,
1488	. u,v,t,q,
1489	. qsurf,
1490	. pcfm, pcfh)
1491	IMPLICIT none
1492	c======================================================================
1493	c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
1494	c (une version strictement identique a l'ancien modele)
1495	c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
1496	c coefficients d'echange turbulent dans l'atmosphere.
1497	c Arguments:
1498	c nsrf-----input-I- indicateur de la nature du sol
1499	c knon-----input-I- nombre de points a traiter
1500	c paprs----input-R- pression a chaque intercouche (en Pa)
1501	c pplay----input-R- pression au milieu de chaque couche (en Pa)
1502	c ts-------input-R- temperature du sol (en Kelvin)
1503	c rugos----input-R- longeur de rugosite (en m)
1504	c u--------input-R- vitesse u
1505	c v--------input-R- vitesse v
1506	c t--------input-R- temperature (K)
1507	c q--------input-R- vapeur d'eau (kg/kg)
1508	c
1509	c itop-----output-I- numero de couche du sommet de la couche limite
1510	c pcfm-----output-R- coefficients a calculer (vitesse)
1511	c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1512	c======================================================================
1513	#include "dimensions.h"
1514	#include "dimphy.h"
1515	#include "YOMCST.h"
1516	#include "indicesol.h"
1517	#include "iniprint.h"
1518	c
1519	c Arguments:
1520	c
1521	INTEGER knon, nsrf
1522	REAL ts(klon)
1523	REAL paprs(klon,klev+1), pplay(klon,klev)
1524	REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
1525	REAL rugos(klon)
1526	c
1527	REAL pcfm(klon,klev), pcfh(klon,klev)
1528	INTEGER itop(klon)
1529	c
1530	c Quelques constantes et options:
1531	c
1532	REAL cepdu2, ckap, cb, cc, cd, clam
1533	PARAMETER (cepdu2 =(0.1)**2)
1534	PARAMETER (CKAP=0.4)
1535	PARAMETER (cb=5.0)
1536	PARAMETER (cc=5.0)
1537	PARAMETER (cd=5.0)
1538	PARAMETER (clam=160.0)
1539	REAL ratqs ! largeur de distribution de vapeur d'eau
1540	PARAMETER (ratqs=0.05)
1541	LOGICAL richum ! utilise le nombre de Richardson humide
1542	PARAMETER (richum=.TRUE.)
1543	REAL ric ! nombre de Richardson critique
1544	PARAMETER(ric=0.4)
1545	REAL prandtl
1546	PARAMETER (prandtl=0.4)
1547	REAL kstable ! diffusion minimale (situation stable)
1548	! GKtest
1549	! PARAMETER (kstable=1.0e-10)
1550	REAL ksta, ksta_ter
1551	cIM: 261103 REAL kstable_ter, kstable_sinon
1552	cIM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6)
1553	cIM: 261103 PARAMETER (kstable_ter = 1.0e-8)
1554	cIM: 261103 PARAMETER (kstable_ter = 1.0e-10)
1555	cIM: 261103 PARAMETER (kstable_sinon = 1.0e-10)
1556	! fin GKtest
1557	REAL mixlen ! constante controlant longueur de melange
1558	PARAMETER (mixlen=35.0)
1559	INTEGER isommet ! le sommet de la couche limite
1560	PARAMETER (isommet=klev)
1561	LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
1562	PARAMETER (tvirtu=.TRUE.)
1563	LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
1564	PARAMETER (opt_ec=.FALSE.)
1565
1566	#include "compbl.h"
1567	c
1568	c Variables locales:
1569	c
1570	INTEGER i, k, kk !IM 120704
1571	REAL zgeop(klon,klev)
1572	REAL zmgeom(klon)
1573	REAL zri(klon)
1574	REAL zl2(klon)
1575
1576	REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
1577	REAL pcfm1(klon), pcfh1(klon)
1578	c
1579	REAL zdphi, zdu2, ztvd, ztvu, zcdn
1580	REAL zscf
1581	REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
1582	REAL z2geomf, zalh2, zalm2, zscfh, zscfm
1583	REAL t_coup
1584	PARAMETER (t_coup=273.15)
1585	cIM
1586	LOGICAL check
1587	PARAMETER (check=.false.)
1588	c
1589	c contre-gradient pour la chaleur sensible: Kelvin/metre
1590	REAL gamt(2:klev)
1591	real qsurf(klon)
1592	c
1593	LOGICAL appel1er
1594	SAVE appel1er
1595	c
1596	c Fonctions thermodynamiques et fonctions d'instabilite
1597	REAL fsta, fins, x
1598	LOGICAL zxli ! utiliser un jeu de fonctions simples
1599	PARAMETER (zxli=.FALSE.)
1600	c
1601	#include "YOETHF.h"
1602	#include "FCTTRE.h"
1603	fsta(x) = 1.0 / (1.0+10.0x(1+8.0*x))
1604	fins(x) = SQRT(1.0-18.0*x)
1605	c
1606	DATA appel1er /.TRUE./
1607	c
1608	IF (appel1er) THEN
1609	if (prt_level > 9) THEN
1610	WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
1611	WRITE(lunout,*)'coefkz, richum:', richum
1612	IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
1613	WRITE(lunout,*)'coefkz, isommet:', isommet
1614	WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
1615	appel1er = .FALSE.
1616	endif
1617	ENDIF
1618	c
1619	c Initialiser les sorties
1620	c
1621	DO k = 1, klev
1622	DO i = 1, knon
1623	pcfm(i,k) = 0.0
1624	pcfh(i,k) = 0.0
1625	ENDDO
1626	ENDDO
1627	DO i = 1, knon
1628	itop(i) = 0
1629	ENDDO
1630
1631	c
1632	c Prescrire la valeur de contre-gradient
1633	c
1634	if (iflag_pbl.eq.1) then
1635	DO k = 3, klev
1636	gamt(k) = -1.0E-03
1637	ENDDO
1638	gamt(2) = -2.5E-03
1639	else
1640	DO k = 2, klev
1641	gamt(k) = 0.0
1642	ENDDO
1643	ENDIF
1644	cIM cf JLD/ GKtest
1645	IF ( nsrf .NE. is_oce ) THEN
1646	cIM 261103 kstable = kstable_ter
1647	kstable = ksta_ter
1648	ELSE
1649	cIM 261103 kstable = kstable_sinon
1650	kstable = ksta
1651	ENDIF
1652	cIM cf JLD/ GKtest fin
1653	c
1654	c Calculer les geopotentiels de chaque couche
1655	c
1656	DO i = 1, knon
1657	zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1658	. * (paprs(i,1)-pplay(i,1))
1659	ENDDO
1660	DO k = 2, klev
1661	DO i = 1, knon
1662	zgeop(i,k) = zgeop(i,k-1)
1663	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1664	. * (pplay(i,k-1)-pplay(i,k))
1665	ENDDO
1666	ENDDO
1667	c
1668	c Calculer le frottement au sol (Cdrag)
1669	c
1670	DO i = 1, knon
1671	u1(i) = u(i,1)
1672	v1(i) = v(i,1)
1673	t1(i) = t(i,1)
1674	q1(i) = q(i,1)
1675	z1(i) = zgeop(i,1)
1676	ENDDO
1677	c
1678	CALL clcdrag(klon, knon, nsrf, zxli,
1679	$ u1, v1, t1, q1, z1,
1680	$ ts, qsurf, rugos,
1681	$ pcfm1, pcfh1)
1682	cIM $ ts, qsurf, rugos,
1683	C
1684	DO i = 1, knon
1685	pcfm(i,1)=pcfm1(i)
1686	pcfh(i,1)=pcfh1(i)
1687	ENDDO
1688	c
1689	c Calculer les coefficients turbulents dans l'atmosphere
1690	c
1691	DO i = 1, knon
1692	itop(i) = isommet
1693	ENDDO
1694
1695
1696	DO k = 2, isommet
1697	DO i = 1, knon
1698	zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
1699	. +(v(i,k)-v(i,k-1))**2)
1700	zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
1701	zdphi =zmgeom(i) / 2.0
1702	zt = (t(i,k)+t(i,k-1)) * 0.5
1703	zq = (q(i,k)+q(i,k-1)) * 0.5
1704	c
1705	c calculer Qs et dQs/dT:
1706	c
1707	IF (thermcep) THEN
1708	zdelta = MAX(0.,SIGN(1.,RTT-zt))
1709	zcvm5 = R5LESRLVTT/RCPD/(1.0+RVTMP2zq)*(1.-zdelta)
1710	. + R5IESRLSTT/RCPD/(1.0+RVTMP2zq)*zdelta
1711	zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
1712	zqs = MIN(0.5,zqs)
1713	zcor = 1./(1.-RETV*zqs)
1714	zqs = zqs*zcor
1715	zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
1716	ELSE
1717	IF (zt .LT. t_coup) THEN
1718	zqs = qsats(zt) / pplay(i,k)
1719	zdqs = dqsats(zt,zqs)
1720	ELSE
1721	zqs = qsatl(zt) / pplay(i,k)
1722	zdqs = dqsatl(zt,zqs)
1723	ENDIF
1724	ENDIF
1725	c
1726	c calculer la fraction nuageuse (processus humide):
1727	c
1728	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
1729	zfr = MAX(0.0,MIN(1.0,zfr))
1730	IF (.NOT.richum) zfr = 0.0
1731	c
1732	c calculer le nombre de Richardson:
1733	c
1734	IF (tvirtu) THEN
1735	ztvd =( t(i,k)
1736	. + zdphi/RCPD/(1.+RVTMP2*zq)
1737	. ( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
1738	. )(1.+RETVq(i,k))
1739	ztvu =( t(i,k-1)
1740	. - zdphi/RCPD/(1.+RVTMP2*zq)
1741	. ( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
1742	. )(1.+RETVq(i,k-1))
1743	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
1744	zri(i) = zri(i)
1745	. + zmgeom(i)zmgeom(i)/RGgamt(k)
1746	. (paprs(i,k)/101325.0)*RKAPPA
1747	. /(zdu20.5(ztvd+ztvu))
1748	c
1749	ELSE ! calcul de Ridchardson compatible LMD5
1750	c
1751	zri(i) =(RCPD*(t(i,k)-t(i,k-1))
1752	. -RD0.5(t(i,k)+t(i,k-1))/paprs(i,k)
1753	. *(pplay(i,k)-pplay(i,k-1))
1754	. )zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
1755	zri(i) = zri(i) +
1756	. zmgeom(i)zmgeom(i)gamt(k)/RG
1757	cSB . /(paprs(i,k)/101325.0)**RKAPPA
1758	. (paprs(i,k)/101325.0)*RKAPPA
1759	. /(zdu20.5(t(i,k-1)+t(i,k)))
1760	ENDIF
1761	c
1762	c finalement, les coefficients d'echange sont obtenus:
1763	c
1764	zcdn=SQRT(zdu2) / zmgeom(i) * RG
1765	c
1766	IF (opt_ec) THEN
1767	z2geomf=zgeop(i,k-1)+zgeop(i,k)
1768	zalm2=(0.5ckap/RGz2geomf
1769	. /(1.+0.5ckap/rg/clamz2geomf))**2
1770	zalh2=(0.5ckap/rgz2geomf
1771	. /(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
1772	IF (zri(i).LT.0.0) THEN ! situation instable
1773	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3
1774	. / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
1775	zscf = SQRT(-zri(i)*zscf)
1776	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
1777	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
1778	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
1779	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
1780	ELSE ! situation stable
1781	zscf=SQRT(1.+cd*zri(i))
1782	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
1783	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
1784	ENDIF
1785	ELSE
1786	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
1787	. /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
1788	pcfm(i,k)=sqrt(max(zcdnzcdn(ric-zri(i))/ric, kstable))
1789	pcfm(i,k)= zl2(i)* pcfm(i,k)
1790	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1791	ENDIF
1792	ENDDO
1793	ENDDO
1794	c
1795	c Au-dela du sommet, pas de diffusion turbulente:
1796	c
1797	DO i = 1, knon
1798	IF (itop(i)+1 .LE. klev) THEN
1799	DO k = itop(i)+1, klev
1800	pcfh(i,k) = 0.0
1801	pcfm(i,k) = 0.0
1802	ENDDO
1803	ENDIF
1804	ENDDO
1805	c
1806	RETURN
1807	END
1808
1809	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t,
1810	. pcfm, pcfh)
1811	IMPLICIT none
1812	c======================================================================
1813	c J'introduit un peu de diffusion sauf dans les endroits
1814	c ou une forte inversion est presente
1815	c On peut dire qu'il represente la convection peu profonde
1816	c
1817	c Arguments:
1818	c nsrf-----input-I- indicateur de la nature du sol
1819	c knon-----input-I- nombre de points a traiter
1820	c paprs----input-R- pression a chaque intercouche (en Pa)
1821	c pplay----input-R- pression au milieu de chaque couche (en Pa)
1822	c t--------input-R- temperature (K)
1823	c
1824	c pcfm-----output-R- coefficients a calculer (vitesse)
1825	c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1826	c======================================================================
1827	#include "dimensions.h"
1828	#include "dimphy.h"
1829	#include "YOMCST.h"
1830	#include "indicesol.h"
1831	#include "iniprint.h"
1832	c
1833	c Arguments:
1834	c
1835	INTEGER knon, nsrf
1836	REAL paprs(klon,klev+1), pplay(klon,klev)
1837	REAL t(klon,klev)
1838	c
1839	REAL pcfm(klon,klev), pcfh(klon,klev)
1840	c
1841	c Quelques constantes et options:
1842	c
1843	REAL prandtl
1844	PARAMETER (prandtl=0.4)
1845	REAL kstable
1846	PARAMETER (kstable=0.002)
1847	ccc PARAMETER (kstable=0.001)
1848	REAL mixlen ! constante controlant longueur de melange
1849	PARAMETER (mixlen=35.0)
1850	REAL seuil ! au-dela l'inversion est consideree trop faible
1851	PARAMETER (seuil=-0.02)
1852	ccc PARAMETER (seuil=-0.04)
1853	ccc PARAMETER (seuil=-0.06)
1854	ccc PARAMETER (seuil=-0.09)
1855	c
1856	c Variables locales:
1857	c
1858	INTEGER i, k, invb(knon)
1859	REAL zl2(knon)
1860	REAL zdthmin(knon), zdthdp
1861	c
1862	c Initialiser les sorties
1863	c
1864	DO k = 1, klev
1865	DO i = 1, knon
1866	pcfm(i,k) = 0.0
1867	pcfh(i,k) = 0.0
1868	ENDDO
1869	ENDDO
1870	c
1871	c Chercher la zone d'inversion forte
1872	c
1873	DO i = 1, knon
1874	invb(i) = klev
1875	zdthmin(i)=0.0
1876	ENDDO
1877	DO k = 2, klev/2-1
1878	DO i = 1, knon
1879	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1880	. - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
1881	zdthdp = zdthdp * 100.0
1882	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1883	. zdthdp.LT.zdthmin(i) ) THEN
1884	zdthmin(i) = zdthdp
1885	invb(i) = k
1886	ENDIF
1887	ENDDO
1888	ENDDO
1889	c
1890	c Introduire une diffusion:
1891	c
1892	DO k = 2, klev
1893	DO i = 1, knon
1894	cIM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean
1895	cIM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
1896	!IM cf JLD/ GKtest TERkz2
1897	! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre
1898	! fin GKtest
1899	IF ( (nsrf.EQ.is_oce) .AND. ! si on est sur ocean et si
1900	. ( (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
1901	. (zdthmin(i).GT.seuil) ) )THEN ! si l'inversion est trop faible
1902	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1))
1903	. /(paprs(i,2)-paprs(i,klev+1)) ))**2
1904	pcfm(i,k)= zl2(i)* kstable
1905	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1906	ENDIF
1907	ENDDO
1908	ENDDO
1909	c
1910	RETURN
1911	END
1912	SUBROUTINE calbeta(dtime,indice,knon,snow,qsol,
1913	. vbeta,vcal,vdif)
1914	IMPLICIT none
1915	c======================================================================
1916	c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
1917	c date: 19940414
1918	c======================================================================
1919	c
1920	c Calculer quelques parametres pour appliquer la couche limite
1921	c ------------------------------------------------------------
1922	#include "dimensions.h"
1923	#include "dimphy.h"
1924	#include "YOMCST.h"
1925	#include "indicesol.h"
1926	#include "iniprint.h"
1927	REAL tau_gl ! temps de relaxation pour la glace de mer
1928	ccc PARAMETER (tau_gl=86400.0*30.0)
1929	PARAMETER (tau_gl=86400.0*5.0)
1930	REAL mx_eau_sol
1931	PARAMETER (mx_eau_sol=150.0)
1932	c
1933	REAL calsol, calsno, calice ! epaisseur du sol: 0.15 m
1934	PARAMETER (calsol=1.0/(2.5578E+06*0.15))
1935	PARAMETER (calsno=1.0/(2.3867E+06*0.15))
1936	PARAMETER (calice=1.0/(5.1444E+06*0.15))
1937	C
1938	INTEGER i
1939	c
1940	REAL dtime
1941	REAL snow(klon), qsol(klon)
1942	INTEGER indice, knon
1943	C
1944	REAL vbeta(klon)
1945	REAL vcal(klon)
1946	REAL vdif(klon)
1947	C
1948
1949	IF (indice.EQ.is_oce) THEN
1950	DO i = 1, knon
1951	vcal(i) = 0.0
1952	vbeta(i) = 1.0
1953	vdif(i) = 0.0
1954	ENDDO
1955	ENDIF
1956	c
1957	IF (indice.EQ.is_sic) THEN
1958	DO i = 1, knon
1959	vcal(i) = calice
1960	IF (snow(i) .GT. 0.0) vcal(i) = calsno
1961	vbeta(i) = 1.0
1962	vdif(i) = 1.0/tau_gl
1963	ccc vdif(i) = calice/tau_gl ! c'etait une erreur
1964	ENDDO
1965	ENDIF
1966	c
1967	IF (indice.EQ.is_ter) THEN
1968	DO i = 1, knon
1969	vcal(i) = calsol
1970	IF (snow(i) .GT. 0.0) vcal(i) = calsno
1971	vbeta(i) = MIN(2.0*qsol(i)/mx_eau_sol, 1.0)
1972	vdif(i) = 0.0
1973	ENDDO
1974	ENDIF
1975	c
1976	IF (indice.EQ.is_lic) THEN
1977	DO i = 1, knon
1978	vcal(i) = calice
1979	IF (snow(i) .GT. 0.0) vcal(i) = calsno
1980	vbeta(i) = 1.0
1981	vdif(i) = 0.0
1982	ENDDO
1983	ENDIF
1984	c
1985	RETURN
1986	END
1987	C======================================================================
1988	SUBROUTINE nonlocal(knon, paprs, pplay,
1989	. tsol,beta,u,v,t,q,
1990	. cd_h, cd_m, pcfh, pcfm, cgh, cgq)
1991	IMPLICIT none
1992	c======================================================================
1993	c Laurent Li (LMD/CNRS), le 30 septembre 1998
1994	c Couche limite non-locale. Adaptation du code du CCM3.
1995	c Code non teste, donc a ne pas utiliser.
1996	c======================================================================
1997	c Nonlocal scheme that determines eddy diffusivities based on a
1998	c diagnosed boundary layer height and a turbulent velocity scale.
1999	c Also countergradient effects for heat and moisture are included.
2000	c
2001	c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
2002	c Local versus nonlocal boundary-layer diffusion in a global climate
2003	c model. J. of Climate, vol. 6, 1825-1842.
2004	c======================================================================
2005	#include "dimensions.h"
2006	#include "dimphy.h"
2007	#include "YOMCST.h"
2008	#include "iniprint.h"
2009	c
2010	c Arguments:
2011	c
2012	INTEGER knon ! nombre de points a calculer
2013	REAL tsol(klon) ! temperature du sol (K)
2014	REAL beta(klon) ! efficacite d'evaporation (entre 0 et 1)
2015	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
2016	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
2017	REAL u(klon,klev) ! vitesse U (m/s)
2018	REAL v(klon,klev) ! vitesse V (m/s)
2019	REAL t(klon,klev) ! temperature (K)
2020	REAL q(klon,klev) ! vapeur d'eau (kg/kg)
2021	REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
2022	REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
2023	c
2024	INTEGER isommet
2025	PARAMETER (isommet=klev)
2026	REAL vk
2027	PARAMETER (vk=0.40)
2028	REAL ricr
2029	PARAMETER (ricr=0.4)
2030	REAL fak
2031	PARAMETER (fak=8.5)
2032	REAL fakn
2033	PARAMETER (fakn=7.2)
2034	REAL onet
2035	PARAMETER (onet=1.0/3.0)
2036	REAL t_coup
2037	PARAMETER(t_coup=273.15)
2038	REAL zkmin
2039	PARAMETER (zkmin=0.01)
2040	REAL betam
2041	PARAMETER (betam=15.0)
2042	REAL betah
2043	PARAMETER (betah=15.0)
2044	REAL betas
2045	PARAMETER (betas=5.0)
2046	REAL sffrac
2047	PARAMETER (sffrac=0.1)
2048	REAL binm
2049	PARAMETER (binm=betam*sffrac)
2050	REAL binh
2051	PARAMETER (binh=betah*sffrac)
2052	REAL ccon
2053	PARAMETER (ccon=faksffracvk)
2054	c
2055	REAL z(klon,klev)
2056	REAL pcfm(klon,klev), pcfh(klon,klev)
2057	c
2058	INTEGER i, k
2059	REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
2060	REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
2061	REAL khfs(klon) ! surface kinematic heat flux [mK/s]
2062	REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
2063	REAL heatv(klon) ! surface virtual heat flux
2064	REAL ustar(klon)
2065	REAL rino(klon,klev) ! bulk Richardon no. from level to ref lev
2066	LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
2067	LOGICAL stblev(klon) ! stable pbl with levels within pbl
2068	LOGICAL unslev(klon) ! unstbl pbl with levels within pbl
2069	LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr
2070	LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr
2071	LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
2072	REAL pblh(klon)
2073	REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
2074	REAL cgq(klon,2:klev) ! counter-gradient term for constituents
2075	REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
2076	REAL obklen(klon)
2077	REAL ztvd, ztvu, zdu2
2078	REAL therm(klon) ! thermal virtual temperature excess
2079	REAL phiminv(klon) ! inverse phi function for momentum
2080	REAL phihinv(klon) ! inverse phi function for heat
2081	REAL wm(klon) ! turbulent velocity scale for momentum
2082	REAL fak1(klon) ! kustarpblh
2083	REAL fak2(klon) ! kwmpblh
2084	REAL fak3(klon) ! fakn*wstr/wm
2085	REAL pblk(klon) ! level eddy diffusivity for momentum
2086	REAL pr(klon) ! Prandtl number for eddy diffusivities
2087	REAL zl(klon) ! zmzp / Obukhov length
2088	REAL zh(klon) ! zmzp / pblh
2089	REAL zzh(klon) ! (1-(zmzp/pblh))**2
2090	REAL wstr(klon) ! w*, convective velocity scale
2091	REAL zm(klon) ! current level height
2092	REAL zp(klon) ! current level height + one level up
2093	REAL zcor, zdelta, zcvm5, zxqs
2094	REAL fac, pblmin, zmzp, term
2095	c
2096	#include "YOETHF.h"
2097	#include "FCTTRE.h"
2098	c
2099	c Initialisation
2100	c
2101	DO i = 1, klon
2102	pcfh(i,1) = cd_h(i)
2103	pcfm(i,1) = cd_m(i)
2104	ENDDO
2105	DO k = 2, klev
2106	DO i = 1, klon
2107	pcfh(i,k) = zkmin
2108	pcfm(i,k) = zkmin
2109	cgs(i,k) = 0.0
2110	cgh(i,k) = 0.0
2111	cgq(i,k) = 0.0
2112	ENDDO
2113	ENDDO
2114	c
2115	c Calculer les hauteurs de chaque couche
2116	c
2117	DO i = 1, knon
2118	z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
2119	. * (paprs(i,1)-pplay(i,1)) / RG
2120	ENDDO
2121	DO k = 2, klev
2122	DO i = 1, knon
2123	z(i,k) = z(i,k-1)
2124	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
2125	. * (pplay(i,k-1)-pplay(i,k)) / RG
2126	ENDDO
2127	ENDDO
2128	c
2129	DO i = 1, knon
2130	IF (thermcep) THEN
2131	zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
2132	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
2133	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
2134	zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
2135	zxqs=MIN(0.5,zxqs)
2136	zcor=1./(1.-retv*zxqs)
2137	zxqs=zxqs*zcor
2138	ELSE
2139	IF (tsol(i).LT.t_coup) THEN
2140	zxqs = qsats(tsol(i)) / paprs(i,1)
2141	ELSE
2142	zxqs = qsatl(tsol(i)) / paprs(i,1)
2143	ENDIF
2144	ENDIF
2145	zx_alf1 = 1.0
2146	zx_alf2 = 1.0 - zx_alf1
2147	zxt = (t(i,1)+z(i,1)RG/RCPD/(1.+RVTMP2q(i,1)))
2148	. (1.+RETVq(i,1))*zx_alf1
2149	. + (t(i,2)+z(i,2)RG/RCPD/(1.+RVTMP2q(i,2)))
2150	. (1.+RETVq(i,2))*zx_alf2
2151	zxu = u(i,1)zx_alf1+u(i,2)zx_alf2
2152	zxv = v(i,1)zx_alf1+v(i,2)zx_alf2
2153	zxq = q(i,1)zx_alf1+q(i,2)zx_alf2
2154	zxmod = 1.0+SQRT(zxu2+zxv2)
2155	khfs(i) = (tsol(i)(1.+RETVq(i,1))-zxt) zxmodcd_h(i)
2156	kqfs(i) = (zxqs-zxq) zxmodcd_h(i) * beta(i)
2157	heatv(i) = khfs(i) + 0.61zxtkqfs(i)
2158	taux = zxu zxmodcd_m(i)
2159	tauy = zxv zxmodcd_m(i)
2160	ustar(i) = SQRT(taux2+tauy2)
2161	ustar(i) = MAX(SQRT(ustar(i)),0.01)
2162	ENDDO
2163	c
2164	DO i = 1, knon
2165	rino(i,1) = 0.0
2166	check(i) = .TRUE.
2167	pblh(i) = z(i,1)
2168	obklen(i) = -t(i,1)ustar(i)3/(RGvk*heatv(i))
2169	ENDDO
2170
2171	C
2172	C PBL height calculation:
2173	C Search for level of pbl. Scan upward until the Richardson number between
2174	C the first level and the current level exceeds the "critical" value.
2175	C
2176	fac = 100.0
2177	DO k = 1, isommet
2178	DO i = 1, knon
2179	IF (check(i)) THEN
2180	zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
2181	zdu2 = max(zdu2,1.0e-20)
2182	ztvd =(t(i,k)+z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,k)))
2183	. (1.+RETVq(i,k))
2184	ztvu =(t(i,1)-z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
2185	. (1.+RETVq(i,1))
2186	rino(i,k) = (z(i,k)-z(i,1))RG(ztvd-ztvu)
2187	. /(zdu20.5(ztvd+ztvu))
2188	IF (rino(i,k).GE.ricr) THEN
2189	pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
2190	. (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
2191	check(i) = .FALSE.
2192	ENDIF
2193	ENDIF
2194	ENDDO
2195	ENDDO
2196
2197	C
2198	C Set pbl height to maximum value where computation exceeds number of
2199	C layers allowed
2200	C
2201	DO i = 1, knon
2202	if (check(i)) pblh(i) = z(i,isommet)
2203	ENDDO
2204	C
2205	C Improve estimate of pbl height for the unstable points.
2206	C Find unstable points (sensible heat flux is upward):
2207	C
2208	DO i = 1, knon
2209	IF (heatv(i) .GT. 0.) THEN
2210	unstbl(i) = .TRUE.
2211	check(i) = .TRUE.
2212	ELSE
2213	unstbl(i) = .FALSE.
2214	check(i) = .FALSE.
2215	ENDIF
2216	ENDDO
2217	C
2218	C For the unstable case, compute velocity scale and the
2219	C convective temperature excess:
2220	C
2221	DO i = 1, knon
2222	IF (check(i)) THEN
2223	phiminv(i) = (1.-binmpblh(i)/obklen(i))*onet
2224	wm(i)= ustar(i)*phiminv(i)
2225	therm(i) = heatv(i)*fak/wm(i)
2226	rino(i,1) = 0.0
2227	ENDIF
2228	ENDDO
2229	C
2230	C Improve pblh estimate for unstable conditions using the
2231	C convective temperature excess:
2232	C
2233	DO k = 1, isommet
2234	DO i = 1, knon
2235	IF (check(i)) THEN
2236	zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
2237	zdu2 = max(zdu2,1.0e-20)
2238	ztvd =(t(i,k)+z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,k)))
2239	. (1.+RETVq(i,k))
2240	ztvu =(t(i,1)+therm(i)-z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
2241	. (1.+RETVq(i,1))
2242	rino(i,k) = (z(i,k)-z(i,1))RG(ztvd-ztvu)
2243	. /(zdu20.5(ztvd+ztvu))
2244	IF (rino(i,k).GE.ricr) THEN
2245	pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
2246	. (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
2247	check(i) = .FALSE.
2248	ENDIF
2249	ENDIF
2250	ENDDO
2251	ENDDO
2252	C
2253	C Set pbl height to maximum value where computation exceeds number of
2254	C layers allowed
2255	C
2256	DO i = 1, knon
2257	if (check(i)) pblh(i) = z(i,isommet)
2258	ENDDO
2259	C
2260	C Points for which pblh exceeds number of pbl layers allowed;
2261	C set to maximum
2262	C
2263	DO i = 1, knon
2264	IF (check(i)) pblh(i) = z(i,isommet)
2265	ENDDO
2266	C
2267	C PBL height must be greater than some minimum mechanical mixing depth
2268	C Several investigators have proposed minimum mechanical mixing depth
2269	C relationships as a function of the local friction velocity, u*. We
2270	C make use of a linear relationship of the form h = c u* where c=700.
2271	C The scaling arguments that give rise to this relationship most often
2272	C represent the coefficient c as some constant over the local coriolis
2273	C parameter. Here we make use of the experimental results of Koracin
2274	C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
2275	C where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
2276	C latitude value for f so that c = 0.07/f = 700.
2277	C
2278	DO i = 1, knon
2279	pblmin = 700.0*ustar(i)
2280	pblh(i) = MAX(pblh(i),pblmin)
2281	ENDDO
2282	C
2283	C pblh is now available; do preparation for diffusivity calculation:
2284	C
2285	DO i = 1, knon
2286	pblk(i) = 0.0
2287	fak1(i) = ustar(i)pblh(i)vk
2288	C
2289	C Do additional preparation for unstable cases only, set temperature
2290	C and moisture perturbations depending on stability.
2291	C
2292	IF (unstbl(i)) THEN
2293	zxt=(t(i,1)-z(i,1)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
2294	. (1.+RETVq(i,1))
2295	phiminv(i) = (1. - binmpblh(i)/obklen(i))*onet
2296	phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
2297	wm(i) = ustar(i)*phiminv(i)
2298	fak2(i) = wm(i)pblh(i)vk
2299	wstr(i) = (heatv(i)RGpblh(i)/zxt)**onet
2300	fak3(i) = fakn*wstr(i)/wm(i)
2301	ENDIF
2302	ENDDO
2303
2304	C Main level loop to compute the diffusivities and
2305	C counter-gradient terms:
2306	C
2307	DO 1000 k = 2, isommet
2308	C
2309	C Find levels within boundary layer:
2310	C
2311	DO i = 1, knon
2312	unslev(i) = .FALSE.
2313	stblev(i) = .FALSE.
2314	zm(i) = z(i,k-1)
2315	zp(i) = z(i,k)
2316	IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
2317	IF (zm(i) .LT. pblh(i)) THEN
2318	zmzp = 0.5*(zm(i) + zp(i))
2319	zh(i) = zmzp/pblh(i)
2320	zl(i) = zmzp/obklen(i)
2321	zzh(i) = 0.
2322	IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
2323	C
2324	C stblev for points zm < plbh and stable and neutral
2325	C unslev for points zm < plbh and unstable
2326	C
2327	IF (unstbl(i)) THEN
2328	unslev(i) = .TRUE.
2329	ELSE
2330	stblev(i) = .TRUE.
2331	ENDIF
2332	ENDIF
2333	ENDDO
2334	C
2335	C Stable and neutral points; set diffusivities; counter-gradient
2336	C terms zero for stable case:
2337	C
2338	DO i = 1, knon
2339	IF (stblev(i)) THEN
2340	IF (zl(i).LE.1.) THEN
2341	pblk(i) = fak1(i)zh(i)zzh(i)/(1. + betas*zl(i))
2342	ELSE
2343	pblk(i) = fak1(i)zh(i)zzh(i)/(betas + zl(i))
2344	ENDIF
2345	pcfm(i,k) = pblk(i)
2346	pcfh(i,k) = pcfm(i,k)
2347	ENDIF
2348	ENDDO
2349	C
2350	C unssrf, unstable within surface layer of pbl
2351	C unsout, unstable within outer layer of pbl
2352	C
2353	DO i = 1, knon
2354	unssrf(i) = .FALSE.
2355	unsout(i) = .FALSE.
2356	IF (unslev(i)) THEN
2357	IF (zh(i).lt.sffrac) THEN
2358	unssrf(i) = .TRUE.
2359	ELSE
2360	unsout(i) = .TRUE.
2361	ENDIF
2362	ENDIF
2363	ENDDO
2364	C
2365	C Unstable for surface layer; counter-gradient terms zero
2366	C
2367	DO i = 1, knon
2368	IF (unssrf(i)) THEN
2369	term = (1. - betamzl(i))*onet
2370	pblk(i) = fak1(i)zh(i)zzh(i)*term
2371	pr(i) = term/sqrt(1. - betah*zl(i))
2372	ENDIF
2373	ENDDO
2374	C
2375	C Unstable for outer layer; counter-gradient terms non-zero:
2376	C
2377	DO i = 1, knon
2378	IF (unsout(i)) THEN
2379	pblk(i) = fak2(i)zh(i)zzh(i)
2380	cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
2381	cgh(i,k) = khfs(i)*cgs(i,k)
2382	pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
2383	cgq(i,k) = kqfs(i)*cgs(i,k)
2384	ENDIF
2385	ENDDO
2386	C
2387	C For all unstable layers, set diffusivities
2388	C
2389	DO i = 1, knon
2390	IF (unslev(i)) THEN
2391	pcfm(i,k) = pblk(i)
2392	pcfh(i,k) = pblk(i)/pr(i)
2393	ENDIF
2394	ENDDO
2395	1000 continue ! end of level loop
2396
2397	RETURN
2398	END

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

Download in other formats:

Original Format