Changeset 2787 for LMDZ5/branches/testing/libf/phylmd/wake.F90

Timestamp:

Jan 30, 2017, 5:54:45 PM (8 years ago)

Author:

Laurent Fairhead

Message:

Merged trunk changes r2727:2785 into testing branch

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phylmd/wake.F90 (modified) (18 diffs)

LMDZ5/branches/testing/libf/phylmd/wake.F90

@@ -2 +2 @@
 ! $Id$
 
-SUBROUTINE wake(p, ph, pi, dtime, &
+SUBROUTINE wake(znatsurf, p, ph, pi, dtime, &
                 te0, qe0, omgb, &
                 dtdwn, dqdwn, amdwn, amup, dta, dqa, &
@@ -124 +124 @@
   ! --------------------
 
+  INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf
   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
   REAL, DIMENSION (klon, klev+1),   INTENT(IN)          :: ph
@@ -167 +168 @@
   LOGICAL, SAVE                                         :: first = .TRUE.
   !$OMP THREADPRIVATE(first)
-  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
+!jyg<
+!!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
+  REAL, SAVE, DIMENSION(2)                              :: wdens_ref
+  REAL, SAVE                                            :: stark, coefgw, alpk
+!>jyg
   REAL, SAVE                                            :: crep_upper, crep_sol  
   !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
+
+  LOGICAL, SAVE                                         :: flag_wk_check_trgl
+  !$OMP THREADPRIVATE(flag_wk_check_trgl)
 
   REAL                                                  :: delta_t_min
@@ -200 +208 @@
   REAL, DIMENSION (klon)                                :: z, dz, hw0
   INTEGER, DIMENSION (klon)                             :: ktop, kupper
+
+  ! Variables liées au test de la forme triangulaire du profil de Delta_theta
+  REAL, DIMENSION (klon)                                :: sum_half_dth
+  REAL, DIMENSION (klon)                                :: dz_half
 
   ! Sub-timestep tendencies and related variables
@@ -312 +324 @@
   alpk=0.25
   CALL getin_p('alpk',alpk)
-  wdens_ref=8.E-12
-  CALL getin_p('wdens_ref',wdens_ref)
+!jyg<
+!!  wdens_ref=8.E-12
+!!  CALL getin_p('wdens_ref',wdens_ref)
+  wdens_ref(1)=8.E-12
+  wdens_ref(2)=8.E-12
+  CALL getin_p('wdens_ref_o',wdens_ref(1))    !wake number per unit area ; ocean
+  CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
+!>jyg
   coefgw=4.
   CALL getin_p('coefgw',coefgw)
@@ -319 +337 @@
   WRITE(*,*) 'stark=', stark
   WRITE(*,*) 'alpk=', alpk
-  WRITE(*,*) 'wdens_ref=', wdens_ref
+!jyg<
+!!  WRITE(*,*) 'wdens_ref=', wdens_ref
+  WRITE(*,*) 'wdens_ref_o=', wdens_ref(1)
+  WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
+!>jyg
   WRITE(*,*) 'coefgw=', coefgw 
+
+  flag_wk_check_trgl=.false.
+  CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl)
+  WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl
 
   first=.false.
@@ -328 +354 @@
   ! Les densites peuvent evoluer si les poches debordent
   ! (voir au tout debut de la boucle sur les substeps)
-  wdens(:) = wdens_ref
+!jyg<
+!!  wdens(:) = wdens_ref
+  DO i = 1,klon
+    wdens(i) = wdens_ref(znatsurf(i)+1)
+  ENDDO 
+!>jyg
 
   ! print*,'stark',stark
@@ -824 +855 @@
 
   IF (prt_level>=10) THEN
-    PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout) ', &
-                      sigmaw(igout), cstar(igout), wape(igout)
+    PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
+                      sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
   ENDIF
 
@@ -846 +877 @@
     END DO
     IF (prt_level>=10) THEN
-      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout) ', &
-                          isubstep,wk_adv(igout),cstar(igout),wape(igout)
+      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
+                          isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
     ENDIF
 
@@ -962 +993 @@
     IF (prt_level>=10) THEN
       PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
-      PRINT *, 'wake-4.2, omgtop(igout) ', omgtop(igout)
+      PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
+                          omgtop(igout), ptop(igout), ktop(igout)
     ENDIF
 
@@ -1609 +1641 @@
 
   IF (prt_level>=10) THEN
-    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout) ', &
-                      sigmaw(igout), cstar(igout), wape(igout)
+    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
+                      sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
   ENDIF
 
@@ -1632 +1664 @@
       sum_thvu(i) = 0.
       sum_dth(i) = 0.
+      sum_half_dth(i) = 0.
       sum_dq(i) = 0.
       sum_rho(i) = 0.
@@ -1646 +1679 @@
       av_dtdwn(i) = 0.
       av_dqdwn(i) = 0.
+
+      dthmin(i) = -delta_t_min
     END IF
   END DO
@@ -1685 +1720 @@
       z(i) = 1.
       dz(i) = 1.
+      dz_half(i) = 1.
       sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
       sum_dth(i) = 0.
@@ -1696 +1732 @@
         ! cc
         dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
+        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*rg)
         IF (dz(i)>0) THEN
           z(i) = z(i) + dz(i)
@@ -1707 +1744 @@
           sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
           sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
+!
+          dthmin(i) = min(dthmin(i), dth(i,k))
+        END IF
+        IF (dz_half(i)>0) THEN
+          sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
         END IF
       END IF
@@ -1750 +1792 @@
 
   ! Filter out bad wakes
+
+  IF (flag_wk_check_trgl) THEN
+    ! Check triangular shape of dth profile
+    DO i = 1, klon
+      IF (ok_qx_qw(i)) THEN
+        !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
+        !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
+        !!                2.*sum_dth(i)/(hw0(i)*dthmin(i))
+        !! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
+        !!                sum_half_dth(i), sum_dth(i)
+        IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
+          wape2(i) = -1.
+          !! print *,'wake, rej 1'
+        ELSE IF (abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
+          wape2(i) = -1.
+          !! print *,'wake, rej 2'
+        ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
+          wape2(i) = -1.
+          !! print *,'wake, rej 3'
+        END IF
+      END IF
+    END DO
+  END IF
+
 
   DO k = 1, klev
@@ -1882 +1948 @@
 
       ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
-      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
+!jyg<
+!!      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
+      IF (ok_qx_qw(i)) THEN
+!>jyg
         ! cc
         dtls(i, k) = dtls(i, k)/dtime