Context Navigation

clmain.F @ 682

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

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

Download in other formats:

Original Format