Changeset 2980 for LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/tropopause_m.F90

Timestamp:

Aug 3, 2017, 10:31:46 AM (7 years ago)

Author:

acozic

Message:

merges with trunk for rev 2968 and 2971 --
2968 = improved method for ozone forcing
2971 = control outputs for debug are removed

File:

• LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/tropopause_m.F90 (modified) (4 diffs)

LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/tropopause_m.F90

@@ -2 +2 @@
 
 !  USE phys_local_var_mod, ONLY: ptrop => pr_tropopause
+  IMPLICIT NONE
   PRIVATE
   PUBLIC :: dyn_tropopause
@@ -9 +10 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE dyn_tropopause(t,ts,paprs,pplay,presnivs,rot,ptrop,tpot0,pvor0,pmin0,pmax0)
+FUNCTION dyn_tropopause(t, ts, paprs, pplay, rot, thet0, pvor0)
 !
 !-------------------------------------------------------------------------------
-  USE assert_m,      ONLY: assert
-  USE assert_eq_m,   ONLY: assert_eq
-  USE comvert_mod,   ONLY: preff
-  USE dimphy,        ONLY: klon, klev
-  USE geometry_mod,  ONLY: latitude_deg, longitude_deg
-  USE interpolation, ONLY: locate
-  IMPLICIT NONE
+  USE assert_m,     ONLY: assert
+  USE assert_eq_m,  ONLY: assert_eq
+  USE dimphy,       ONLY: klon, klev
+  USE geometry_mod, ONLY: latitude_deg, longitude_deg
+  USE vertical_layers_mod, ONLY: aps, bps, preff
+
 !-------------------------------------------------------------------------------
 ! Arguments:
+  REAL ::     dyn_tropopause(klon) !--- Pressure at tropopause
   REAL, INTENT(IN)  ::      t(:,:) !--- Cells-centers temperature
   REAL, INTENT(IN)  ::     ts(:)   !--- Surface       temperature
   REAL, INTENT(IN)  ::  paprs(:,:) !--- Cells-edges   pressure
   REAL, INTENT(IN)  ::  pplay(:,:) !--- Cells-centers pressure
-  REAL, INTENT(IN)  :: presnivs(:) !--- Cells-centers nominal pressure
   REAL, INTENT(IN)  ::    rot(:,:) !--- Cells-centers relative vorticity
-  REAL, INTENT(OUT) :: ptrop(klon) !--- Tropopause pressure
-  REAL, INTENT(IN), OPTIONAL :: tpot0, pvor0, pmin0, pmax0
+  REAL, INTENT(IN), OPTIONAL :: thet0, pvor0
 !-------------------------------------------------------------------------------
 ! Local variables:
   include "YOMCST.h"
   CHARACTER(LEN=80)  :: sub
-  INTEGER :: it, i, nsrf, k, nh, kmin, kmax
-  REAL    :: p1, p2, tpo0, pvo0, pmn0, pmx0, al
-  REAL, DIMENSION(klon,klev)   :: t_edg, t_pot
-  REAL, DIMENSION(klon,klev-1) :: pvort
+  INTEGER :: i, k, kb, kt, kp, ib, ie, nw
+  REAL    :: al, th0, pv0
+  REAL,    DIMENSION(klon,klev) :: tpot_cen, tpot_edg, pvor_cen
+  REAL,    PARAMETER :: sg0=0.75  !--- Start level for PV=cte search loop
+  INTEGER, PARAMETER :: nadj=3    !--- Adjacent levs nb for thresholds detection
+  REAL,    PARAMETER :: w(5)=[0.1,0.25,0.3,0.25,0.1] !--- Vertical smoothing
+  INTEGER, SAVE :: k0
+  LOGICAL, SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(k0,first)
 !-------------------------------------------------------------------------------
   sub='dyn_tropopause'
@@ -42 +46 @@
   CALL assert(SIZE(t ,2)==klev, TRIM(sub)//" t klev")
   CALL assert(SIZE(ts,1)==klon, TRIM(sub)//" ts klon")
-  CALL assert(SIZE(presnivs,1)==klev, TRIM(sub)//" presnivs klev")
   CALL assert(SHAPE(paprs)==[klon,klev+1],TRIM(sub)//" paprs shape")
   CALL assert(SHAPE(pplay)==[klon,klev  ],TRIM(sub)//" pplay shape")
@@ -48 +51 @@
 
   !--- DEFAULT THRESHOLDS
-  tpo0=380.;   IF(PRESENT(tpot0)) tpo0=tpot0  !--- In kelvins
-  pvo0=2.0;    IF(PRESENT(pvor0)) pvo0=pvor0  !--- In PVU
-  pmn0= 8000.; IF(PRESENT(pmin0)) pmn0=pmin0  !--- In pascals
-  pmx0=50000.; IF(PRESENT(pmax0)) pmx0=pmax0  !--- In pascals
-  kmin=klev-locate(presnivs(klev:1:-1),pmx0)+1
-  kmax=klev-locate(presnivs(klev:1:-1),pmn0)+1
+  th0=380.; IF(PRESENT(thet0)) th0=thet0   !--- In kelvins
+  pv0=  2.; IF(PRESENT(pvor0)) pv0=pvor0   !--- In PVU
+  IF(first) THEN
+    DO k0=1,klev; IF(aps(k0)/preff+bps(k0)<sg0) EXIT; END DO; first=.FALSE.
+  END IF
 
-  !--- POTENTIAL TEMPERATURE AT CELLS CENTERS
-  DO k= 1, klev
-    DO i = 1, klon
-      t_pot(i,k) = t(i,k)*(preff/pplay(i,k))**RKAPPA
+  !--- POTENTIAL TEMPERATURE AT CELLS CENTERS AND INTERFACES
+  DO i = 1,klon
+    tpot_cen(i,1) = t(i,1)*(preff/pplay(i,1))**RKAPPA
+    tpot_edg(i,1) = ts(i) *(preff/paprs(i,1))**RKAPPA
+    DO k=2,klev
+      al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
+      tpot_cen(i,k) =  t(i,k)                        *(preff/pplay(i,k))**RKAPPA
+      tpot_edg(i,k) = (t(i,k-1)+al*(t(i,k)-t(i,k-1)))*(preff/paprs(i,k))**RKAPPA
+      !--- FORCE QUANTITIES TO BE GROWING
+      IF(tpot_edg(i,k)<tpot_edg(i,k-1)) tpot_edg(i,k)=tpot_edg(i,k-1)+1.E-5
+      IF(tpot_cen(i,k)<tpot_cen(i,k-1)) tpot_cen(i,k)=tpot_cen(i,k-1)+1.E-5
     END DO
-  END DO
-
-  !--- TEMPERATURE AT CELLS INTERFACES (except in top layer)
-  t_edg(:,1) = ts(:)
-  DO k= 2, klev
-    DO i = 1, klon
-      al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
-      t_edg(i,k) = t(i,k) + al*(t(i,k-1)-t(i,k))
-    END DO
+    !--- VERTICAL SMOOTHING
+    tpot_cen(i,:)=smooth(tpot_cen(i,:),w)
+    tpot_edg(i,:)=smooth(tpot_edg(i,:),w)
   END DO
 
   !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer)
-  DO k= 1, klev-1
-    DO i = 1, klon
-      IF(paprs(i,k)==paprs(i,k+1)) THEN; pvort(i,k)=HUGE(1.); CYCLE; END IF
-      pvort(i,k) = -1.E6 * RG * 2.*ROMEGA*ABS(SIN(latitude_deg(i)*RPI/180.))   &
-                         * ( t_edg(i,k+1)*(preff/paprs(i,k+1) )**RKAPPA        &
-                           - t_edg(i,k  )*(preff/paprs(i,k  ) )**RKAPPA)       &
-                         / ( paprs(i,k+1)  -     paprs(i,k  ) )
+  DO i = 1, klon
+    DO k= 1, klev-1
+      pvor_cen(i,k)=-1.E6*RG*(rot(i,k)+2.*ROMEGA*SIN(latitude_deg(i)*RPI/180.))&
+                   * (tpot_edg(i,k+1)-tpot_edg(i,k)) / (paprs(i,k+1)-paprs(i,k))
     END DO
+    !--- VERTICAL SMOOTHING
+    pvor_cen(i,1:klev-1)=smooth(pvor_cen(i,1:klev-1),w)
   END DO
 
-  !--- LOCATE TROPOPAUSE
-  ptrop(:)=0.
+  !--- LOCATE TROPOPAUSE: LOWEST POINT BETWEEN THETA=380K AND PV=2PVU SURFACES.
   DO i = 1, klon
-    !--- Dynamical tropopause
-    it=kmax; DO WHILE(pvort(i,it)>pvo0.AND.it>=kmin); it=it-1; END DO
-    IF(pvort(i,it+1)/=pvort(i,it).AND.ALL([kmax,kmin-1]/=it) &
-      .AND.ALL([pvort(i,it),pvort(i,it+1)]/=HUGE(1.))) THEN
-      al = (pvort(i,it+1)-pvo0)/(pvort(i,it+1)-pvort(i,it))
-      ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
-    END IF
-    !--- Potential temperature iso-surface
-    it = kmin-1+locate(t_pot(i,kmin:kmax),tpo0)
-    IF(t_pot(i,it+1)/=t_pot(i,it).AND.ALL([kmin-1,kmax]/=it)) THEN
-      al = (t_pot(i,it+1)-tpo0)/(t_pot(i,it+1)-t_pot(i,it))
-      ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
-    END IF
+    !--- UPPER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM TOP
+    DO kt=klev-1,1,-1; IF(ALL(ABS(pvor_cen(i,kt-nadj:kt))<=pv0)) EXIT; END DO
+    !--- LOWER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM BOTTOM
+    DO kb=k0,klev-1;   IF(ALL(ABS(pvor_cen(i,kb:kb+nadj))> pv0)) EXIT; END DO; kb=kb-1
+    !--- ISO-THETA POINT: THETA=380K       STARTING FROM TOP
+    DO kp=klev-1,1,-1; IF(ALL(ABS(tpot_cen(i,kp-nadj:kp))<=th0)) EXIT; END DO
+    !--- CHOOSE BETWEEN LOWER AND UPPER TROPOPAUSE
+    IF(2*COUNT(ABS(pvor_cen(i,kb:kt))>pv0)>kt-kb+1) kt=kb
+    !--- PV-DEFINED TROPOPAUSE
+    al = (ABS(pvor_cen(i,kt+1))-pv0)/ABS(pvor_cen(i,kt+1)-pvor_cen(i,kt))
+    dyn_tropopause(i) = pplay(i,kt+1)*(pplay(i,kt)/pplay(i,kt+1))**al
+    !--- THETA=380K IN THE TROPICAL REGION
+    al = (tpot_cen(i,kp+1)-th0)/(tpot_cen(i,kp+1)-tpot_cen(i,kp))
+    dyn_tropopause(i) = MAX( pplay(i,kp+1)*(pplay(i,kp)/pplay(i,kp+1))**al,    &
+                            dyn_tropopause(i) )
   END DO
 
-END SUBROUTINE dyn_tropopause
+END FUNCTION dyn_tropopause
+
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION smooth(v,w)
+!
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)         :: v(:), w(:)
+  REAL, DIMENSION(SIZE(v)) :: smooth
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: nv, nw, k, kb, ke, lb, le
+!-------------------------------------------------------------------------------
+  nv=SIZE(v); nw=(SIZE(w)-1)/2
+  DO k=1,nv
+    kb=MAX(k-nw,1 ); lb=MAX(2+nw   -k,1)
+    ke=MIN(k+nw,nv); le=MIN(1+nw+nv-k,1+2*nw)
+    smooth(k)=SUM(v(kb:ke)*w(lb:le))/SUM(w(lb:le))
+  END DO
+
+END FUNCTION smooth
+!
+!-------------------------------------------------------------------------------
 
 END MODULE tropopause_m