source: LMDZ6/trunk/libf/phylmd/lmdz_wake_pkupper.f90

Last change on this file was 5808, checked in by fhourdin, 2 months ago

Sommet des poches froides
File size: 11.5 KB
Line 
1MODULE lmdz_wake_pkupper
2PUBLIC wake_pkupper
3CONTAINS
4
5SUBROUTINE wake_pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
6                    dth, hw_, rho, delta_t_min_in, &
7                    ktop, wk_adv, h_zzz, ptop1, ktop1)
8
9USE lmdz_wake_ini , ONLY : wk_pupper
10USE lmdz_wake_ini , ONLY : RG
11USE lmdz_wake_ini , ONLY : hwmin
12USE lmdz_wake_ini , ONLY : iflag_wk_new_ptop, wk_delta_t_min, wk_frac_int_delta_t
13USE lmdz_wake_ini , ONLY : wk_int_delta_t_min
14
15IMPLICIT NONE
16
17INTEGER,                              INTENT(IN) :: klon,klev
18REAL,       DIMENSION (klon,klev+1) , INTENT(IN) :: ph, p
19REAL,       DIMENSION (klon,klev+1) , INTENT(IN) :: rho
20LOGICAL,    DIMENSION (klon)        , INTENT(IN) :: wk_adv
21REAL,       DIMENSION (klon,klev+1) , INTENT(IN) :: dth
22REAL,                                 INTENT(IN) :: delta_t_min_in
23
24
25REAL,       DIMENSION (klon)  , INTENT(OUT)        :: hw_
26REAL,       DIMENSION (klon)  , INTENT(OUT)        :: ptop
27INTEGER,    DIMENSION (klon)  , INTENT(OUT)        :: Ktop
28REAL,       DIMENSION (klon)  , INTENT(OUT)        :: pupper
29INTEGER,    DIMENSION (klon)  , INTENT(OUT)        :: kupper
30REAL,       DIMENSION (klon)  , INTENT(OUT)        :: h_zzz       !!
31REAL,       DIMENSION (klon)  , INTENT(OUT)        :: Ptop1      !!
32INTEGER,    DIMENSION (klon)  , INTENT(OUT)        :: ktop1      !!
33
34INTEGER :: i,k
35
36LOGICAL,    DIMENSION (klon)       :: wk_active
37REAL                               :: delta_t_min
38REAL,     DIMENSION (klon)         :: dthmin
39REAL,     DIMENSION (klon)         :: ptop_provis,ptop_new
40REAL,     DIMENSION (klon)         :: z, dz
41REAL,     DIMENSION (klon)         :: sum_dth
42
43INTEGER,     DIMENSION (klon)                     :: k_ptop_provis
44REAL,     DIMENSION (klon)                     :: zk_ptop_provis
45REAL,     DIMENSION (klon)                        :: omega        !!
46REAL,     DIMENSION (klon,klev+1)                 :: int_dth      !!
47REAL,     DIMENSION (klon,klev+1)                 :: dzz          !!
48REAL,     DIMENSION (klon,klev+1)                 :: zzz          !!
49REAL,     DIMENSION (klon)                 :: frac_int_dth          !!
50REAL                                              :: ddd!!
51
52REAL :: www
53
54INTEGER, SAVE :: ipas=0
55!$OMP THREADPRIVATE(ipas)
56
57
58
59!INTEGER, SAVE :: compte=0
60
61! LJYF : a priori z, dz sum_dth sont aussi des variables internes
62! Les eliminer apres verification convergence numerique
63
64!compte=compte+1
65!print*,'compte=',compte
66
67    ! Determine Ptop from buoyancy integral
68    ! ---------------------------------------
69
70    ! -     1/ Pressure of the level where dth changes sign.
71    !print*,'WAKE LJYF'
72
73
74if (iflag_wk_new_ptop==0) then
75    delta_t_min=delta_t_min_in
76else
77    delta_t_min=wk_delta_t_min
78endif
79
80    DO i = 1, klon
81        ptop_provis(i) = ph(i, 1)
82        k_ptop_provis(i) = 1
83    END DO
84
85    DO k = 2, klev
86      DO i = 1, klon
87        IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
88! LJYF changer :           dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
89            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
90          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
91                            (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
92          k_ptop_provis(i) = k
93        END IF
94      END DO
95    END DO
96
97
98
99    ! -     2/ dth integral
100
101    DO i = 1, klon
102      IF (wk_adv(i)) THEN !!! nrlmd
103        sum_dth(i) = 0.
104        dthmin(i) = -delta_t_min
105        z(i) = 0.
106      END IF
107    END DO
108
109    DO k = 1, klev
110      DO i = 1, klon
111        IF (wk_adv(i)) THEN
112          dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG)
113          IF (dz(i)>0) THEN
114            z(i) = z(i) + dz(i)
115            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
116            dthmin(i) = amin1(dthmin(i), dth(i,k))
117          END IF
118        END IF
119      END DO
120    END DO
121
122    ! -     3/ height of triangle with area= sum_dth and base = dthmin
123
124    DO i = 1, klon
125      IF (wk_adv(i)) THEN
126        hw_(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
127        hw_(i) = amax1(hwmin, hw_(i))
128      END IF
129    END DO
130
131    ! -     4/ now, get Ptop
132
133    DO i = 1, klon
134      IF (wk_adv(i)) THEN !!! nrlmd
135        ktop(i) = 0
136        z(i) = 0.
137      END IF
138    END DO
139
140    DO k = 1, klev
141      DO i = 1, klon
142        IF (wk_adv(i)) THEN
143          dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw_(i)-z(i))
144          IF (dz(i)>0) THEN
145            z(i) = z(i) + dz(i)
146            ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i)
147            ktop(i) = k
148          END IF
149        END IF
150      END DO
151    END DO
152
153    ! 4.5/Correct ktop and ptop
154
155    DO i = 1, klon
156        ptop_new(i) = ptop(i)
157    END DO
158
159    DO k = klev, 2, -1
160      DO i = 1, klon
161        ! IM v3JYG; IF (k .GE. ktop(i)
162        IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
163! LJYF changer :           dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
164            dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
165          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
166                         (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
167        END IF
168      END DO
169    END DO
170
171
172    DO i = 1, klon
173        ptop(i) = ptop_new(i)
174    END DO
175
176    DO k = klev, 1, -1
177      DO i = 1, klon
178        IF (wk_adv(i)) THEN !!! nrlmd
179          IF (ph(i,k+1)<ptop(i)) ktop(i) = k
180        END IF
181      END DO
182    END DO
183 
184!  IF (prt_level>=10) THEN
185!    PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
186!  ENDIF
187
188    ! -----------------------------------------------------------------------
189    ! nouveau calcul de hw et ptop
190    ! -----------------------------------------------------------------------
191!if (iflag_wk_new_ptop>0) then
192do i=1,klon
193   ptop1(i)=ph(i,1)
194   ktop1(i)=1
195   h_zzz(i)=0.
196enddo
197   
198IF (iflag_wk_new_ptop/=0) THEN
199   
200    int_dth(1:klon,1:klev+1)=0.
201    DO i = 1, klon
202       IF (wk_adv(i)) THEN
203          int_dth(i,1) = 0.
204      END IF
205    END DO
206   
207    if (abs(iflag_wk_new_ptop) == 1 ) then
208        DO k = 2, klev+1
209           Do i = 1, klon
210              IF (wk_adv(i)) THEN
211                 if (k<=k_ptop_provis(i)) then
212                      ddd=dth(i,k-1)*(ph(i,k-1) - max(ptop_provis(i),ph(i,k)))
213                      !ddd=dth(i,k-1)*(ph(i,k-1) - ph(i,k))
214                 else
215                      ddd=0.
216                 endif             
217                 int_dth(i,k) = int_dth(i,k-1) + ddd
218              !ELSE
219              !   int_dth(i,k) = 0.
220              END IF
221           END DO
222        END DO
223    else
224        k_ptop_provis(:)=klev+1
225        dthmin(:)=dth(:,1)
226        ! calcul de l'int??grale de dT * dP jusqu'au dernier
227        ! niveau avec dT<0. (en s'assurant qu'on a bien un
228        ! dT negatif plus bas)
229        DO k = 1, klev
230           DO i = 1, klon
231              dthmin(i)=min(dthmin(i),dth(i,k))
232              ddd=dth(i,k)*(ph(i,k)-ph(i,k+1))
233              if (dthmin(i)<0.) then
234                  if (k>=k_ptop_provis(i)) then
235                      ddd=0.
236                  else if (dth(i,k)>=0.) then
237                      ddd=0.
238                      k_ptop_provis(i)=k+1
239                  endif
240              endif
241              int_dth(i,k+1) = int_dth(i,k)+ ddd
242           ENDDO
243        ENDDO
244
245        DO i = 1, klon
246           if ( k_ptop_provis(i)==klev+1 .or. .not. wk_adv(i)) then
247                k_ptop_provis(i)=1
248           endif
249        ENDDO
250    endif ! (abs(iflag_wk_new_ptop) == 1 )
251   ! print*, 'xxx, int_dth', (k,int_dth(1,k),k=1,klev)
252   ! print*, 'xxx, k_ptop_provis', k_ptop_provis(1)
253   
254
255 
256    ! On se limite ?? des poches avec integrale dT * dp < -wk_int_delta_t_min
257    do i=1,klon
258          if (int_dth(i,k_ptop_provis(i)) > -wk_int_delta_t_min .or. k_ptop_provis(i)==1) then
259          !if (1==0) then
260             wk_active(i)=.false.
261             ptop(i)=ph(i,1)
262             ktop(i)=1
263             hw_(i)=0.
264          else
265             wk_active(i)=wk_adv(i)
266          endif
267    enddo
268
269    DO i=1,klon
270       IF (wk_active(i)) THEN
271        frac_int_dth(i)=wk_frac_int_delta_t*int_dth(i,k_ptop_provis(i))
272       ENDIF
273    ENDDO
274    DO k = 1,klev
275       DO i =1, klon
276          IF (wk_active(i)) THEN
277            IF (int_dth(i,k)>=frac_int_dth(i)) THEN
278              ktop1(i) = min(k, k_ptop_provis(i))
279              !ktop1(i) = k
280              !print*,ipas,'yyy ktop1= ',ktop1
281            ENDIF
282          ENDIF
283       END DO
284    END DO
285    !print*, 'LAMINE'
286   
287    DO i = 1, klon
288       IF (wk_active(i)) THEN
289           !print*, ipas,'xxx1, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ',ktop1
290           ddd=int_dth(i,ktop1(i)+1)-int_dth(i,ktop1(i))
291           if (ddd==0.) then
292              omega(i)=0.
293           else
294              omega(i) = (frac_int_dth(i) - int_dth(i,ktop1(i)))/ddd
295           endif
296           !! print*,'OMEGA ',omega(i)
297       END IF
298    END DO
299   
300    !! print*, 'xxx'
301    DO i = 1, klon
302       IF (wk_active(i)) THEN
303      ! print*, 'xxx, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ', &
304      !               int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1)
305      ! print*, 'xxx, omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) ',  &
306      !e               omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1)
307          ptop1(i) = min((1 - omega(i))*ph(i,ktop1(i)) + omega(i)*ph(i,ktop1(i)+1), ph(i,1))
308      END IF
309    END DO
310   
311    DO i=1, klon
312       IF (wk_active(i)) THEN
313           zzz(i, 1) = 0
314       END IF
315     END DO
316    DO k = 1, klev
317       DO i = 1, klon
318           IF (wk_active(i)) THEN         
319              dzz(i,k) = (ph(i,k) - ph(i,k+1))/(rho(i,k)*RG)
320              zzz(i,k+1) = zzz(i,k) + dzz(i,k)
321           END IF
322       END DO
323    END DO
324   
325    DO i =1, klon
326       IF (wk_active(i)) THEN
327           h_zzz(i) = max((1- omega(i))*zzz(i,ktop1(i)) + omega(i)*zzz(i,ktop1(i)+1), hwmin)
328       END IF
329    END DO
330
331
332ENDIF ! (iflag_wk_new_ptop/=0)
333
334!if (iflag_wk_new_ptop==2) then
335IF (iflag_wk_new_ptop>0) THEN
336   do i=1,klon
337      ptop(i)=ptop1(i)
338      ktop(i)=ktop1(i)
339      hw_(i)=h_zzz(i)
340   enddo
341
342!endif
343ENDIF
344
345 kupper = 0
346 
347IF (0.<wk_pupper .and. wk_pupper<1.) THEN
348 ! Choose an integration bound well above wake top
349  ! -----------------------------------------------------------------
350
351  ! Pupper = 50000.  ! melting level
352  ! Pupper = 60000.
353  ! Pupper = 80000.  ! essais pour case_e
354  DO i = 1, klon
355  !  pupper(i) = 0.6*ph(i, 1)
356    pupper(i) = wk_pupper*ph(i, 1)
357    pupper(i) = max(pupper(i), 45000.)
358    ! cc        Pupper(i) = 60000.
359  END DO
360
361ELSE IF (1.<=wk_pupper) THEN
362  DO i=1, klon
363     ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1)
364     !  pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-50.)
365      pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-5000.)
366  END DO
367ELSE
368  www=abs(wk_pupper)
369  DO i=1, klon
370     ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1)
371     !  pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-50.)
372      pupper(i) = min(   exp(www*log(ptop(i))+(1.-www)*log(ph(i, 1)) ) , ptop(i)-5000.)
373  END DO
374END IF
375 
376  ! -5/ Determination de kupper
377
378  DO k = klev, 1, -1
379    DO i = 1, klon
380      IF (ph(i,k+1)<pupper(i)) kupper(i) = k
381    END DO
382  END DO
383
384  ! On evite kupper = 1 et kupper = klev
385  DO i = 1, klon
386    kupper(i) = max(kupper(i), 2)
387    kupper(i) = min(kupper(i), klev-1)
388  END DO
389  !---------- FIN nouveau calcul hw et ptop -------------------------------------
390
391IF (iflag_wk_new_ptop==999) THEN
392    DO i = 1, klon
393    hw_(i)=0.
394    ptop(i)=ph(i,1)
395    Ktop(i)=1
396    pupper(i)=ph(i,2)
397    kupper(i)=2
398    h_zzz(i)=0.
399    Ptop1(i)=ph(i,1)
400    ENDDO
401ENDIF
402
403zk_ptop_provis=k_ptop_provis
404
405    RETURN
406END SUBROUTINE wake_pkupper
407END MODULE lmdz_wake_pkupper
Note: See TracBrowser for help on using the repository browser.