source: LMDZ6/branches/Amaury_dev/libf/phylmd/stdlevvar_mod.F90

Last change on this file was 5160, checked in by abarral, 3 months ago

Put .h into modules
  • Property svn:executable set to *
File size: 25.8 KB
Line 
1MODULE stdlevvar_mod
2
3  ! This module contains main procedures for calculation
4  ! of temperature, specific humidity and wind at a reference level
5
6  USE cdrag_mod
7  USE lmdz_screenp, ONLY: screenp
8  USE screenc_mod
9  IMPLICIT NONE
10
11CONTAINS
12
13  !****************************************************************************************
14
15  !r original routine svn3623
16
17  SUBROUTINE stdlevvar(klon, knon, nsrf, zxli, &
18          u1, v1, t1, q1, z1, &
19          ts1, qsurf, z0m, z0h, psol, pat1, &
20          t_2m, q_2m, t_10m, q_10m, u_10m, ustar, s_pblh, prain, tsol)
21    USE lmdz_flux_arp, ONLY: fsens, flat, betaevap, ust, tg, ok_flux_surf, ok_prescr_ust, ok_prescr_beta, ok_forc_tsurf
22    USE lmdz_yoethf
23    USE lmdz_yomcst
24
25    IMPLICIT NONE
26    !-------------------------------------------------------------------------
27
28    ! Objet : calcul de la temperature et l'humidite relative a 2m et du
29    !         module du vent a 10m a partir des relations de Dyer-Businger et
30    !         des equations de Louis.
31
32    ! Reference : Hess, Colman et McAvaney (1995)
33
34    ! I. Musat, 01.07.2002
35
36    !AM On rajoute en sortie t et q a 10m pr le calcule d'hbtm2 dans clmain
37
38    !-------------------------------------------------------------------------
39
40    ! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
41    ! knon----input-I- nombre de points pour un type de surface
42    ! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
43    ! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
44    ! u1------input-R- vent zonal au 1er niveau du modele
45    ! v1------input-R- vent meridien au 1er niveau du modele
46    ! t1------input-R- temperature de l'air au 1er niveau du modele
47    ! q1------input-R- humidite relative au 1er niveau du modele
48    ! z1------input-R- geopotentiel au 1er niveau du modele
49    ! ts1-----input-R- temperature de l'air a la surface
50    ! qsurf---input-R- humidite relative a la surface
51    ! z0m, z0h---input-R- rugosite
52    ! psol----input-R- pression au sol
53    ! pat1----input-R- pression au 1er niveau du modele
54
55    ! t_2m---output-R- temperature de l'air a 2m
56    ! q_2m---output-R- humidite relative a 2m
57    ! u_10m--output-R- vitesse du vent a 10m
58    !AM
59    ! t_10m--output-R- temperature de l'air a 10m
60    ! q_10m--output-R- humidite specifique a 10m
61    ! ustar--output-R- u*
62
63    INTEGER, INTENT(IN) :: klon, knon, nsrf
64    LOGICAL, INTENT(IN) :: zxli
65    REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, t1, q1, z1, ts1
66    REAL, DIMENSION(klon), INTENT(IN) :: qsurf
67    REAL, DIMENSION(klon), INTENT(INOUT) :: z0m, z0h
68    REAL, DIMENSION(klon), INTENT(IN) :: psol, pat1
69
70    REAL, DIMENSION(klon), INTENT(OUT) :: t_2m, q_2m, ustar
71    REAL, DIMENSION(klon), INTENT(OUT) :: u_10m, t_10m, q_10m
72    REAL, DIMENSION(klon), INTENT(INOUT) :: s_pblh
73    REAL, DIMENSION(klon), INTENT(IN) :: prain
74    REAL, DIMENSION(klon), INTENT(IN) :: tsol
75    !-------------------------------------------------------------------------
76
77    ! Quelques constantes et options:
78
79    ! RKAR : constante de von Karman
80    REAL, PARAMETER :: RKAR = 0.40
81    ! niter : nombre iterations calcul "corrector"
82    !     INTEGER, parameter :: niter=6, ncon=niter-1
83    INTEGER, parameter :: niter = 2, ncon = niter - 1
84
85    ! Variables locales
86    INTEGER :: i, n
87    REAL :: zref
88    REAL, DIMENSION(klon) :: speed
89    ! tpot : temperature potentielle
90    REAL, DIMENSION(klon) :: tpot
91    REAL, DIMENSION(klon) :: zri1, cdran
92    REAL, DIMENSION(klon) :: cdram, cdrah
93    ! ri1 : nb. de Richardson entre la surface --> la 1ere couche
94    REAL, DIMENSION(klon) :: ri1
95    REAL, DIMENSION(klon) :: testar, qstar
96    REAL, DIMENSION(klon) :: zdte, zdq
97    ! lmon : longueur de Monin-Obukhov selon Hess, Colman and McAvaney
98    DOUBLE PRECISION, DIMENSION(klon) :: lmon
99    DOUBLE PRECISION, parameter :: eps = 1.0D-20
100    REAL, DIMENSION(klon) :: delu, delte, delq
101    REAL, DIMENSION(klon) :: u_zref, te_zref, q_zref
102    REAL, DIMENSION(klon) :: temp, pref
103    LOGICAL :: okri
104    REAL, DIMENSION(klon) :: u_zref_p, te_zref_p, temp_p, q_zref_p
105    !convertgence
106    REAL, DIMENSION(klon) :: te_zref_con, q_zref_con
107    REAL, DIMENSION(klon) :: u_zref_c, te_zref_c, temp_c, q_zref_c
108    REAL, DIMENSION(klon) :: ok_pred, ok_corr, zri_zero
109    !     REAL, DIMENSION(klon) :: conv_te, conv_q
110    !-------------------------------------------------------------------------
111    DO i = 1, knon
112      speed(i) = SQRT(u1(i)**2 + v1(i)**2)
113      ri1(i) = 0.0
114    ENDDO
115
116    okri = .FALSE.
117    !      CALL coefcdrag(klon, knon, nsrf, zxli, &
118    ! &                   speed, t1, q1, z1, psol, &
119    ! &                   ts1, qsurf, rugos, okri, ri1,  &
120    ! &                   cdram, cdrah, cdran, zri1, pref)
121    ! Fuxing WANG, 04/03/2015, replace the coefcdrag by the merged version: cdrag
122
123    CALL cdrag(knon, nsrf, &
124            speed, t1, q1, z1, &
125            psol, s_pblh, ts1, qsurf, z0m, z0h, &
126            zri_zero, 0, &
127            cdram, cdrah, zri1, pref, prain, tsol, pat1)
128
129    ! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
130    IF (ok_prescr_ust) THEN
131      DO i = 1, knon
132        PRINT *, 'cdram avant=', cdram(i)
133        cdram(i) = ust * ust / speed(i) / speed(i)
134        PRINT *, 'cdram ust speed apres=', cdram(i), ust, speed
135      ENDDO
136    ENDIF
137
138    !---------Star variables----------------------------------------------------
139
140    DO i = 1, knon
141      ri1(i) = zri1(i)
142      tpot(i) = t1(i) * (psol(i) / pat1(i))**RKAPPA
143      ustar(i) = sqrt(cdram(i) * speed(i) * speed(i))
144      zdte(i) = tpot(i) - ts1(i)
145      zdq(i) = max(q1(i), 0.0) - max(qsurf(i), 0.0)
146
147
148      !IM BUG BUG BUG       zdte(i) = max(zdte(i),1.e-10)
149      zdte(i) = sign(max(abs(zdte(i)), 1.e-10), zdte(i))
150
151      testar(i) = (cdrah(i) * zdte(i) * speed(i)) / ustar(i)
152      qstar(i) = (cdrah(i) * zdq(i) * speed(i)) / ustar(i)
153      lmon(i) = (ustar(i) * ustar(i) * tpot(i)) / &
154              (RKAR * RG * testar(i))
155    ENDDO
156
157    !----------First aproximation of variables at zref --------------------------
158    zref = 2.0
159    CALL screenp(klon, knon, speed, tpot, q1, &
160            ts1, qsurf, z0m, lmon, &
161            ustar, testar, qstar, zref, &
162            delu, delte, delq)
163
164    DO i = 1, knon
165      u_zref(i) = delu(i)
166      q_zref(i) = max(qsurf(i), 0.0) + delq(i)
167      te_zref(i) = ts1(i) + delte(i)
168      temp(i) = te_zref(i) * (psol(i) / pat1(i))**(-RKAPPA)
169      q_zref_p(i) = q_zref(i)
170      !       te_zref_p(i) = te_zref(i)
171      temp_p(i) = temp(i)
172    ENDDO
173
174    ! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
175
176    DO n = 1, niter
177
178      okri = .TRUE.
179      CALL screenc(klon, knon, nsrf, zxli, &
180              u_zref, temp, q_zref, zref, &
181              ts1, qsurf, z0m, z0h, psol, &
182              ustar, testar, qstar, okri, ri1, &
183              pref, delu, delte, delq, s_pblh, prain, tsol, pat1)
184
185      DO i = 1, knon
186        u_zref(i) = delu(i)
187        q_zref(i) = delq(i) + max(qsurf(i), 0.0)
188        te_zref(i) = delte(i) + ts1(i)
189
190        ! return to normal temperature
191
192        temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
193        !         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
194        !                 (1 + RVTMP2 * max(q_zref(i),0.0))
195
196        !IM +++
197        !         IF(temp(i).GT.350.) THEN
198        !           WRITE(*,*) 'temp(i) GT 350 K !!',i,nsrf,temp(i)
199        !         ENDIF
200        !IM ---
201
202        IF(n==ncon) THEN
203          te_zref_con(i) = te_zref(i)
204          q_zref_con(i) = q_zref(i)
205        ENDIF
206
207      ENDDO
208
209    ENDDO
210
211    ! verifier le critere de convergence : 0.25% pour te_zref et 5% pour qe_zref
212
213    !       DO i = 1, knon
214    !         conv_te(i) = (te_zref(i) - te_zref_con(i))/te_zref_con(i)
215    !         conv_q(i) = (q_zref(i) - q_zref_con(i))/q_zref_con(i)
216    !IM +++
217    !         IF(abs(conv_te(i)).GE.0.0025.AND.abs(conv_q(i)).GE.0.05) THEN
218    !           PRINT*,'DIV','i=',i,te_zref_con(i),te_zref(i),conv_te(i), &
219    !           q_zref_con(i),q_zref(i),conv_q(i)
220    !         ENDIF
221    !IM ---
222    !       ENDDO
223
224    DO i = 1, knon
225      q_zref_c(i) = q_zref(i)
226      temp_c(i) = temp(i)
227
228      !       IF(zri1(i).LT.0.) THEN
229      !         IF(nsrf.EQ.1) THEN
230      !           ok_pred(i)=1.
231      !           ok_corr(i)=0.
232      !         ELSE
233      !           ok_pred(i)=0.
234      !           ok_corr(i)=1.
235      !         ENDIF
236      !       ELSE
237      !         ok_pred(i)=0.
238      !         ok_corr(i)=1.
239      !       ENDIF
240
241      ok_pred(i) = 0.
242      ok_corr(i) = 1.
243
244      t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
245      q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
246      !IM +++
247      !       IF(n.EQ.niter) THEN
248      !       IF(t_2m(i).LT.t1(i).AND.t_2m(i).LT.ts1(i)) THEN
249      !         PRINT*,' BAD t2m LT ',i,nsrf,t_2m(i),t1(i),ts1(i)
250      !       ELSEIF(t_2m(i).GT.t1(i).AND.t_2m(i).GT.ts1(i)) THEN
251      !         PRINT*,' BAD t2m GT ',i,nsrf,t_2m(i),t1(i),ts1(i)
252      !       ENDIF
253      !       ENDIF
254      !IM ---
255    ENDDO
256
257
258    !----------First aproximation of variables at zref --------------------------
259
260    zref = 10.0
261    CALL screenp(klon, knon, speed, tpot, q1, ts1, qsurf, z0m, lmon, ustar, &
262            testar, qstar, zref, delu, delte, delq)
263
264    DO i = 1, knon
265      u_zref(i) = delu(i)
266      q_zref(i) = max(qsurf(i), 0.0) + delq(i)
267      te_zref(i) = ts1(i) + delte(i)
268      temp(i) = te_zref(i) * (psol(i) / pat1(i))**(-RKAPPA)
269      !       temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
270      !                 (1 + RVTMP2 * max(q_zref(i),0.0))
271      u_zref_p(i) = u_zref(i)
272    ENDDO
273
274    ! Iteration of the variables at the reference level zref : corrector ; see Hess & McAvaney, 1995
275
276    DO n = 1, niter
277
278      okri = .TRUE.
279      CALL screenc(klon, knon, nsrf, zxli, &
280              u_zref, temp, q_zref, zref, &
281              ts1, qsurf, z0m, z0h, psol, &
282              ustar, testar, qstar, okri, ri1, &
283              pref, delu, delte, delq, s_pblh, prain, tsol, pat1)
284
285      DO i = 1, knon
286        u_zref(i) = delu(i)
287        q_zref(i) = delq(i) + max(qsurf(i), 0.0)
288        te_zref(i) = delte(i) + ts1(i)
289        temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
290        !         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
291        !                   (1 + RVTMP2 * max(q_zref(i),0.0))
292      ENDDO
293
294    ENDDO
295
296    DO i = 1, knon
297      u_zref_c(i) = u_zref(i)
298
299      u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
300
301      !AM
302      q_zref_c(i) = q_zref(i)
303      temp_c(i) = temp(i)
304      t_10m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
305      q_10m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
306      !MA
307    ENDDO
308
309  END SUBROUTINE stdlevvar
310
311  SUBROUTINE stdlevvarn(klon, knon, nsrf, zxli, &
312          u1, v1, t1, q1, z1, &
313          ts1, qsurf, z0m, z0h, psol, pat1, &
314          t_2m, q_2m, t_10m, q_10m, u_10m, ustar, &
315          n2mout)
316
317    USE lmdz_ioipsl_getin_p, ONLY: getin_p
318    USE lmdz_flux_arp, ONLY: fsens, flat, betaevap, ust, tg, ok_flux_surf, ok_prescr_ust, ok_prescr_beta, ok_forc_tsurf
319    USE lmdz_yoethf
320    USE lmdz_yomcst
321
322    IMPLICIT NONE
323    !-------------------------------------------------------------------------
324
325    ! Objet : calcul de la temperature et l'humidite relative a 2m et du
326    !         module du vent a 10m a partir des relations de Dyer-Businger et
327    !         des equations de Louis.
328
329    ! Reference : Hess, Colman et McAvaney (1995)
330
331    ! I. Musat, 01.07.2002
332
333    !AM On rajoute en sortie t et q a 10m pr le calcule d'hbtm2 dans clmain
334
335    !-------------------------------------------------------------------------
336
337    ! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
338    ! knon----input-I- nombre de points pour un type de surface
339    ! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
340    ! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
341    ! u1------input-R- vent zonal au 1er niveau du modele
342    ! v1------input-R- vent meridien au 1er niveau du modele
343    ! t1------input-R- temperature de l'air au 1er niveau du modele
344    ! q1------input-R- humidite relative au 1er niveau du modele
345    ! z1------input-R- geopotentiel au 1er niveau du modele
346    ! ts1-----input-R- temperature de l'air a la surface
347    ! qsurf---input-R- humidite relative a la surface
348    ! z0m, z0h---input-R- rugosite
349    ! psol----input-R- pression au sol
350    ! pat1----input-R- pression au 1er niveau du modele
351
352    ! t_2m---output-R- temperature de l'air a 2m
353    ! q_2m---output-R- humidite relative a 2m
354    ! u_2m--output-R- vitesse du vent a 2m
355    ! u_10m--output-R- vitesse du vent a 10m
356    ! ustar--output-R- u*
357    !AM
358    ! t_10m--output-R- temperature de l'air a 10m
359    ! q_10m--output-R- humidite specifique a 10m
360
361    INTEGER, INTENT(IN) :: klon, knon, nsrf
362    LOGICAL, INTENT(IN) :: zxli
363    REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, t1, q1, ts1, z1
364    REAL, DIMENSION(klon), INTENT(INOUT) :: z0m, z0h
365    REAL, DIMENSION(klon), INTENT(IN) :: qsurf
366    REAL, DIMENSION(klon), INTENT(IN) :: psol, pat1
367
368    REAL, DIMENSION(klon), INTENT(OUT) :: t_2m, q_2m, ustar
369    REAL, DIMENSION(klon), INTENT(OUT) :: u_10m, t_10m, q_10m
370    INTEGER, DIMENSION(klon, 6), INTENT(OUT) :: n2mout
371
372    REAL, DIMENSION(klon) :: u_2m
373    REAL, DIMENSION(klon) :: cdrm2m, cdrh2m, ri2m
374    REAL, DIMENSION(klon) :: cdram, cdrah, zri1
375    REAL, DIMENSION(klon) :: cdmn1, cdhn1, fm1, fh1
376    REAL, DIMENSION(klon) :: cdmn2m, cdhn2m, fm2m, fh2m
377    REAL, DIMENSION(klon) :: ri2m_new
378    REAL, DIMENSION(klon) :: s_pblh
379    REAL, DIMENSION(klon) :: prain
380    REAL, DIMENSION(klon) :: tsol
381    !-------------------------------------------------------------------------
382
383    ! Quelques constantes et options:
384
385    ! RKAR : constante de von Karman
386    REAL, PARAMETER :: RKAR = 0.40
387    ! niter : nombre iterations calcul "corrector"
388    !     INTEGER, parameter :: niter=6, ncon=niter-1
389    !IM 071020     INTEGER, parameter :: niter=2, ncon=niter-1
390    INTEGER, parameter :: niter = 2, ncon = niter
391    !     INTEGER, parameter :: niter=6, ncon=niter
392
393    ! Variables locales
394    INTEGER :: i, n
395    REAL :: zref
396    REAL, DIMENSION(klon) :: speed
397    ! tpot : temperature potentielle
398    REAL, DIMENSION(klon) :: tpot
399    REAL, DIMENSION(klon) :: cdran
400    ! ri1 : nb. de Richardson entre la surface --> la 1ere couche
401    REAL, DIMENSION(klon) :: ri1
402    DOUBLE PRECISION, parameter :: eps = 1.0D-20
403    REAL, DIMENSION(klon) :: delu, delte, delq
404    REAL, DIMENSION(klon) :: delh, delm
405    REAL, DIMENSION(klon) :: delh_new, delm_new
406    REAL, DIMENSION(klon) :: u_zref, te_zref, q_zref
407    REAL, DIMENSION(klon) :: u_zref_pnew, te_zref_pnew, q_zref_pnew
408    REAL, DIMENSION(klon) :: temp, pref
409    REAL, DIMENSION(klon) :: temp_new, pref_new
410    LOGICAL :: okri
411    REAL, DIMENSION(klon) :: u_zref_p, te_zref_p, temp_p, q_zref_p
412    REAL, DIMENSION(klon) :: u_zref_p_new, te_zref_p_new, temp_p_new, q_zref_p_new
413    !convergence
414    REAL, DIMENSION(klon) :: te_zref_con, q_zref_con
415    REAL, DIMENSION(klon) :: u_zref_c, te_zref_c, temp_c, q_zref_c
416    REAL, DIMENSION(klon) :: ok_pred, ok_corr
417
418    REAL, DIMENSION(klon) :: cdrm10m, cdrh10m, ri10m
419    REAL, DIMENSION(klon) :: cdmn10m, cdhn10m, fm10m, fh10m
420    REAL, DIMENSION(klon) :: cdn2m, cdn1, zri_zero
421    REAL :: CEPDUE, zdu2
422    INTEGER :: nzref, nz1
423    LOGICAL, DIMENSION(klon) :: ok_t2m_toosmall, ok_t2m_toobig
424    LOGICAL, DIMENSION(klon) :: ok_q2m_toosmall, ok_q2m_toobig
425    LOGICAL, DIMENSION(klon) :: ok_u2m_toobig
426    LOGICAL, DIMENSION(klon) :: ok_t10m_toosmall, ok_t10m_toobig
427    LOGICAL, DIMENSION(klon) :: ok_q10m_toosmall, ok_q10m_toobig
428    LOGICAL, DIMENSION(klon) :: ok_u10m_toobig
429    INTEGER, DIMENSION(klon, 6) :: n10mout
430
431    !-------------------------------------------------------------------------
432    CEPDUE = 0.1
433
434    ! n2mout : compteur des pas de temps ou t2m,q2m ou u2m sont en dehors des intervalles
435    !          [tsurf, temp], [qsurf, q1] ou [0, speed]
436    ! n10mout : compteur des pas de temps ou t10m,q10m ou u10m sont en dehors des intervalles
437    !          [tsurf, temp], [qsurf, q1] ou [0, speed]
438
439    n2mout(:, :) = 0
440    n10mout(:, :) = 0
441
442    DO i = 1, knon
443      speed(i) = MAX(SQRT(u1(i)**2 + v1(i)**2), CEPDUE)
444      ri1(i) = 0.0
445    ENDDO
446
447    okri = .FALSE.
448    CALL cdrag(knon, nsrf, &
449            speed, t1, q1, z1, &
450            psol, s_pblh, ts1, qsurf, z0m, z0h, &
451            zri_zero, 0, &
452            cdram, cdrah, zri1, pref, prain, tsol, pat1)
453
454    DO i = 1, knon
455      ri1(i) = zri1(i)
456      tpot(i) = t1(i) * (psol(i) / pat1(i))**RKAPPA
457      zdu2 = MAX(CEPDUE * CEPDUE, speed(i)**2)
458      ustar(i) = sqrt(cdram(i) * zdu2)
459
460    ENDDO
461
462    !----------First aproximation of variables at zref --------------------------
463    zref = 2.0
464
465    ! calcul first-guess en utilisant dans les calculs à 2m
466    ! le Richardson de la premiere couche atmospherique
467
468    CALL screencn(klon, knon, nsrf, zxli, &
469            speed, tpot, q1, zref, &
470            ts1, qsurf, z0m, z0h, psol, &
471            cdram, cdrah, okri, &
472            ri1, 1, &
473            pref_new, delm_new, delh_new, ri2m, &
474            s_pblh, prain, tsol, pat1)
475
476    DO i = 1, knon
477      u_zref(i) = delm_new(i) * speed(i)
478      u_zref_p(i) = u_zref(i)
479      q_zref(i) = delh_new(i) * max(q1(i), 0.0) + &
480              max(qsurf(i), 0.0) * (1 - delh_new(i))
481      q_zref_p(i) = q_zref(i)
482      te_zref(i) = delh_new(i) * tpot(i) + ts1(i) * (1 - delh_new(i))
483      te_zref_p(i) = te_zref(i)
484
485      ! return to normal temperature
486      temp(i) = te_zref(i) * (psol(i) / pref_new(i))**(-RKAPPA)
487      temp_p(i) = temp(i)
488
489      ! compteurs ici
490
491      ok_t2m_toosmall(i) = te_zref(i)<tpot(i).AND. &
492              te_zref(i)<ts1(i)
493      ok_t2m_toobig(i) = te_zref(i)>tpot(i).AND. &
494              te_zref(i)>ts1(i)
495      ok_q2m_toosmall(i) = q_zref(i)<q1(i).AND. &
496              q_zref(i)<qsurf(i)
497      ok_q2m_toobig(i) = q_zref(i)>q1(i).AND. &
498              q_zref(i)>qsurf(i)
499      ok_u2m_toobig(i) = u_zref(i)>speed(i)
500
501      IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i)) THEN
502        n2mout(i, 1) = n2mout(i, 1) + 1
503      ENDIF
504      IF(ok_q2m_toosmall(i).OR.ok_q2m_toobig(i)) THEN
505        n2mout(i, 3) = n2mout(i, 3) + 1
506      ENDIF
507      IF(ok_u2m_toobig(i)) THEN
508        n2mout(i, 5) = n2mout(i, 5) + 1
509      ENDIF
510
511      IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i).OR. &
512              ok_q2m_toosmall(i).OR.ok_q2m_toobig(i).OR. &
513              ok_u2m_toobig(i)) THEN
514        delm_new(i) = min(max(delm_new(i), 0.), 1.)
515        delh_new(i) = min(max(delh_new(i), 0.), 1.)
516        u_zref(i) = delm_new(i) * speed(i)
517        u_zref_p(i) = u_zref(i)
518        q_zref(i) = delh_new(i) * max(q1(i), 0.0) + &
519                max(qsurf(i), 0.0) * (1 - delh_new(i))
520        q_zref_p(i) = q_zref(i)
521        te_zref(i) = delh_new(i) * tpot(i) + ts1(i) * (1 - delh_new(i))
522        te_zref_p(i) = te_zref(i)
523
524        ! return to normal temperature
525        temp(i) = te_zref(i) * (psol(i) / pref_new(i))**(-RKAPPA)
526        temp_p(i) = temp(i)
527      ENDIF
528
529    ENDDO
530
531    ! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
532
533    DO n = 1, niter
534
535      okri = .TRUE.
536      CALL screencn(klon, knon, nsrf, zxli, &
537              u_zref, temp, q_zref, zref, &
538              ts1, qsurf, z0m, z0h, psol, &
539              cdram, cdrah, okri, &
540              ri1, 0, &
541              pref, delm, delh, ri2m, &
542              s_pblh, prain, tsol, pat1)
543
544      DO i = 1, knon
545        u_zref(i) = delm(i) * speed(i)
546        q_zref(i) = delh(i) * max(q1(i), 0.0) + &
547                max(qsurf(i), 0.0) * (1 - delh(i))
548        te_zref(i) = delh(i) * tpot(i) + ts1(i) * (1 - delh(i))
549
550        ! return to normal temperature
551        temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
552
553        ! compteurs ici
554
555        ok_t2m_toosmall(i) = te_zref(i)<tpot(i).AND. &
556                te_zref(i)<ts1(i)
557        ok_t2m_toobig(i) = te_zref(i)>tpot(i).AND. &
558                te_zref(i)>ts1(i)
559        ok_q2m_toosmall(i) = q_zref(i)<q1(i).AND. &
560                q_zref(i)<qsurf(i)
561        ok_q2m_toobig(i) = q_zref(i)>q1(i).AND. &
562                q_zref(i)>qsurf(i)
563        ok_u2m_toobig(i) = u_zref(i)>speed(i)
564
565        IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i)) THEN
566          n2mout(i, 2) = n2mout(i, 2) + 1
567        ENDIF
568        IF(ok_q2m_toosmall(i).OR.ok_q2m_toobig(i)) THEN
569          n2mout(i, 4) = n2mout(i, 4) + 1
570        ENDIF
571        IF(ok_u2m_toobig(i)) THEN
572          n2mout(i, 6) = n2mout(i, 6) + 1
573        ENDIF
574
575        IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i).OR. &
576                ok_q2m_toosmall(i).OR.ok_q2m_toobig(i).OR. &
577                ok_u2m_toobig(i)) THEN
578          delm(i) = min(max(delm(i), 0.), 1.)
579          delh(i) = min(max(delh(i), 0.), 1.)
580          u_zref(i) = delm(i) * speed(i)
581          q_zref(i) = delh(i) * max(q1(i), 0.0) + &
582                  max(qsurf(i), 0.0) * (1 - delh(i))
583          te_zref(i) = delh(i) * tpot(i) + ts1(i) * (1 - delh(i))
584          temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
585        ENDIF
586
587        IF(n==ncon) THEN
588          te_zref_con(i) = te_zref(i)
589          q_zref_con(i) = q_zref(i)
590        ENDIF
591
592      ENDDO
593
594    ENDDO
595
596    DO i = 1, knon
597      q_zref_c(i) = q_zref(i)
598      temp_c(i) = temp(i)
599
600      ok_pred(i) = 0.
601      ok_corr(i) = 1.
602
603      t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
604      q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
605
606      u_zref_c(i) = u_zref(i)
607      u_2m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
608    ENDDO
609
610
611    !----------First aproximation of variables at zref --------------------------
612
613    zref = 10.0
614
615    CALL screencn(klon, knon, nsrf, zxli, &
616            speed, tpot, q1, zref, &
617            ts1, qsurf, z0m, z0h, psol, &
618            cdram, cdrah, okri, &
619            ri1, 1, &
620            pref_new, delm_new, delh_new, ri10m, &
621            s_pblh, prain, tsol, pat1)
622
623    DO i = 1, knon
624      u_zref(i) = delm_new(i) * speed(i)
625      q_zref(i) = delh_new(i) * max(q1(i), 0.0) + &
626              max(qsurf(i), 0.0) * (1 - delh_new(i))
627      te_zref(i) = delh_new(i) * tpot(i) + ts1(i) * (1 - delh_new(i))
628      temp(i) = te_zref(i) * (psol(i) / pref_new(i))**(-RKAPPA)
629      u_zref_p(i) = u_zref(i)
630
631      ! compteurs ici
632
633      ok_t10m_toosmall(i) = te_zref(i)<tpot(i).AND. &
634              te_zref(i)<ts1(i)
635      ok_t10m_toobig(i) = te_zref(i)>tpot(i).AND. &
636              te_zref(i)>ts1(i)
637      ok_q10m_toosmall(i) = q_zref(i)<q1(i).AND. &
638              q_zref(i)<qsurf(i)
639      ok_q10m_toobig(i) = q_zref(i)>q1(i).AND. &
640              q_zref(i)>qsurf(i)
641      ok_u10m_toobig(i) = u_zref(i)>speed(i)
642
643      IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i)) THEN
644        n10mout(i, 1) = n10mout(i, 1) + 1
645      ENDIF
646      IF(ok_q10m_toosmall(i).OR.ok_q10m_toobig(i)) THEN
647        n10mout(i, 3) = n10mout(i, 3) + 1
648      ENDIF
649      IF(ok_u10m_toobig(i)) THEN
650        n10mout(i, 5) = n10mout(i, 5) + 1
651      ENDIF
652
653      IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i).OR. &
654              ok_q10m_toosmall(i).OR.ok_q10m_toobig(i).OR. &
655              ok_u10m_toobig(i)) THEN
656        delm_new(i) = min(max(delm_new(i), 0.), 1.)
657        delh_new(i) = min(max(delh_new(i), 0.), 1.)
658        u_zref(i) = delm_new(i) * speed(i)
659        u_zref_p(i) = u_zref(i)
660        q_zref(i) = delh_new(i) * max(q1(i), 0.0) + &
661                max(qsurf(i), 0.0) * (1 - delh_new(i))
662        te_zref(i) = delh_new(i) * tpot(i) + ts1(i) * (1 - delh_new(i))
663        temp(i) = te_zref(i) * (psol(i) / pref_new(i))**(-RKAPPA)
664      ENDIF
665
666    ENDDO
667
668    ! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
669
670    DO n = 1, niter
671
672      okri = .TRUE.
673      CALL screencn(klon, knon, nsrf, zxli, &
674              u_zref, temp, q_zref, zref, &
675              ts1, qsurf, z0m, z0h, psol, &
676              cdram, cdrah, okri, &
677              ri1, 0, &
678              pref, delm, delh, ri10m, &
679              s_pblh, prain, tsol, pat1)
680
681      DO i = 1, knon
682        u_zref(i) = delm(i) * speed(i)
683        q_zref(i) = delh(i) * max(q1(i), 0.0) + &
684                max(qsurf(i), 0.0) * (1 - delh(i))
685        te_zref(i) = delh(i) * tpot(i) + ts1(i) * (1 - delh(i))
686
687        ! return to normal temperature
688        temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
689
690        ! compteurs ici
691
692        ok_t10m_toosmall(i) = te_zref(i)<tpot(i).AND. &
693                te_zref(i)<ts1(i)
694        ok_t10m_toobig(i) = te_zref(i)>tpot(i).AND. &
695                te_zref(i)>ts1(i)
696        ok_q10m_toosmall(i) = q_zref(i)<q1(i).AND. &
697                q_zref(i)<qsurf(i)
698        ok_q10m_toobig(i) = q_zref(i)>q1(i).AND. &
699                q_zref(i)>qsurf(i)
700        ok_u10m_toobig(i) = u_zref(i)>speed(i)
701
702        IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i)) THEN
703          n10mout(i, 2) = n10mout(i, 2) + 1
704        ENDIF
705        IF(ok_q10m_toosmall(i).OR.ok_q10m_toobig(i)) THEN
706          n10mout(i, 4) = n10mout(i, 4) + 1
707        ENDIF
708        IF(ok_u10m_toobig(i)) THEN
709          n10mout(i, 6) = n10mout(i, 6) + 1
710        ENDIF
711
712        IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i).OR. &
713                ok_q10m_toosmall(i).OR.ok_q10m_toobig(i).OR. &
714                ok_u10m_toobig(i)) THEN
715          delm(i) = min(max(delm(i), 0.), 1.)
716          delh(i) = min(max(delh(i), 0.), 1.)
717          u_zref(i) = delm(i) * speed(i)
718          q_zref(i) = delh(i) * max(q1(i), 0.0) + &
719                  max(qsurf(i), 0.0) * (1 - delh(i))
720          te_zref(i) = delh(i) * tpot(i) + ts1(i) * (1 - delh(i))
721          temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
722        ENDIF
723
724        IF(n==ncon) THEN
725          te_zref_con(i) = te_zref(i)
726          q_zref_con(i) = q_zref(i)
727        ENDIF
728
729      ENDDO
730
731    ENDDO
732
733    DO i = 1, knon
734      q_zref_c(i) = q_zref(i)
735      temp_c(i) = temp(i)
736
737      ok_pred(i) = 0.
738      ok_corr(i) = 1.
739
740      t_10m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
741      q_10m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
742
743      u_zref_c(i) = u_zref(i)
744      u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
745    ENDDO
746
747  END SUBROUTINE stdlevvarn
748
749END MODULE stdlevvar_mod
Note: See TracBrowser for help on using the repository browser.