Context Navigation

clmain.F @ 3919

Last change on this file since 3919 was 704, checked in by Laurent Fairhead, 18 years ago
Inclusion des modifs de Y. Meurdesoif pour la version V3 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 76.7 KB

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

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

Download in other formats:

Original Format