source: LMDZ5/trunk/libf/phylmd/tropopause_m.F90 @ 3605

Last change on this file since 3605 was 2992, checked in by oboucher, 7 years ago

Adding p, z, and t of tropopause in output files.
Stratosphere_mask is streamlined a bit.
This includes a treatment of missing values.

File size: 5.3 KB
Line 
1MODULE tropopause_m
2
3  IMPLICIT NONE
4  PRIVATE
5  PUBLIC :: dyn_tropopause
6
7CONTAINS
8
9!-------------------------------------------------------------------------------
10!
11FUNCTION dyn_tropopause(t, ts, paprs, pplay, rot, thet0, pvor0)
12!
13!-------------------------------------------------------------------------------
14  USE assert_m,     ONLY: assert
15  USE assert_eq_m,  ONLY: assert_eq
16  USE dimphy,       ONLY: klon, klev
17  USE geometry_mod, ONLY: latitude_deg, longitude_deg
18  USE vertical_layers_mod, ONLY: aps, bps, preff
19
20!-------------------------------------------------------------------------------
21! Arguments:
22  REAL ::     dyn_tropopause(klon) !--- Pressure at tropopause
23  REAL, INTENT(IN)  ::      t(:,:) !--- Cells-centers temperature
24  REAL, INTENT(IN)  ::     ts(:)   !--- Surface       temperature
25  REAL, INTENT(IN)  ::  paprs(:,:) !--- Cells-edges   pressure
26  REAL, INTENT(IN)  ::  pplay(:,:) !--- Cells-centers pressure
27  REAL, INTENT(IN)  ::    rot(:,:) !--- Cells-centers relative vorticity
28  REAL, INTENT(IN), OPTIONAL :: thet0, pvor0
29!-------------------------------------------------------------------------------
30! Local variables:
31  include "YOMCST.h"
32  CHARACTER(LEN=80)  :: sub
33  INTEGER :: i, k, kb, kt, kp, ib, ie, nw
34  REAL    :: al, th0, pv0
35  REAL,    DIMENSION(klon,klev) :: tpot_cen, tpot_edg, pvor_cen
36  REAL,    PARAMETER :: sg0=0.75  !--- Start level for PV=cte search loop
37  INTEGER, PARAMETER :: nadj=3    !--- Adjacent levs nb for thresholds detection
38  REAL,    PARAMETER :: w(5)=[0.1,0.25,0.3,0.25,0.1] !--- Vertical smoothing
39  INTEGER, SAVE :: k0
40  LOGICAL, SAVE :: first=.TRUE.
41!$OMP THREADPRIVATE(k0,first)
42!-------------------------------------------------------------------------------
43  sub='dyn_tropopause'
44  CALL assert(SIZE(t ,1)==klon, TRIM(sub)//" t klon")
45  CALL assert(SIZE(t ,2)==klev, TRIM(sub)//" t klev")
46  CALL assert(SIZE(ts,1)==klon, TRIM(sub)//" ts klon")
47  CALL assert(SHAPE(paprs)==[klon,klev+1],TRIM(sub)//" paprs shape")
48  CALL assert(SHAPE(pplay)==[klon,klev  ],TRIM(sub)//" pplay shape")
49  CALL assert(SHAPE(rot)  ==[klon,klev  ],TRIM(sub)//" rot shape")
50
51  !--- DEFAULT THRESHOLDS
52  th0=380.; IF(PRESENT(thet0)) th0=thet0   !--- In kelvins
53  pv0=  2.; IF(PRESENT(pvor0)) pv0=pvor0   !--- In PVU
54  IF(first) THEN
55    DO k0=1,klev; IF(aps(k0)/preff+bps(k0)<sg0) EXIT; END DO; first=.FALSE.
56  END IF
57
58  !--- POTENTIAL TEMPERATURE AT CELLS CENTERS AND INTERFACES
59  DO i = 1,klon
60    tpot_cen(i,1) = t(i,1)*(preff/pplay(i,1))**RKAPPA
61    tpot_edg(i,1) = ts(i) *(preff/paprs(i,1))**RKAPPA
62    DO k=2,klev
63      al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
64      tpot_cen(i,k) =  t(i,k)                        *(preff/pplay(i,k))**RKAPPA
65      tpot_edg(i,k) = (t(i,k-1)+al*(t(i,k)-t(i,k-1)))*(preff/paprs(i,k))**RKAPPA
66      !--- FORCE QUANTITIES TO BE GROWING
67      IF(tpot_edg(i,k)<tpot_edg(i,k-1)) tpot_edg(i,k)=tpot_edg(i,k-1)+1.E-5
68      IF(tpot_cen(i,k)<tpot_cen(i,k-1)) tpot_cen(i,k)=tpot_cen(i,k-1)+1.E-5
69    END DO
70    !--- VERTICAL SMOOTHING
71    tpot_cen(i,:)=smooth(tpot_cen(i,:),w)
72    tpot_edg(i,:)=smooth(tpot_edg(i,:),w)
73  END DO
74
75  !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer)
76  DO i = 1, klon
77    DO k= 1, klev-1
78      pvor_cen(i,k)=-1.E6*RG*(rot(i,k)+2.*ROMEGA*SIN(latitude_deg(i)*RPI/180.))&
79                   * (tpot_edg(i,k+1)-tpot_edg(i,k)) / (paprs(i,k+1)-paprs(i,k))
80    END DO
81    !--- VERTICAL SMOOTHING
82    pvor_cen(i,1:klev-1)=smooth(pvor_cen(i,1:klev-1),w)
83  END DO
84
85  !--- LOCATE TROPOPAUSE: LOWEST POINT BETWEEN THETA=380K AND PV=2PVU SURFACES.
86  DO i = 1, klon
87    !--- UPPER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM TOP
88    DO kt=klev-1,1,-1; IF(ALL(ABS(pvor_cen(i,kt-nadj:kt))<=pv0)) EXIT; END DO
89    !--- LOWER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM BOTTOM
90    DO kb=k0,klev-1;   IF(ALL(ABS(pvor_cen(i,kb:kb+nadj))> pv0)) EXIT; END DO; kb=kb-1
91    !--- ISO-THETA POINT: THETA=380K       STARTING FROM TOP
92    DO kp=klev-1,1,-1; IF(ALL(ABS(tpot_cen(i,kp-nadj:kp))<=th0)) EXIT; END DO
93    !--- CHOOSE BETWEEN LOWER AND UPPER TROPOPAUSE
94    IF(2*COUNT(ABS(pvor_cen(i,kb:kt))>pv0)>kt-kb+1) kt=kb
95    !--- PV-DEFINED TROPOPAUSE
96    al = (ABS(pvor_cen(i,kt+1))-pv0)/ABS(pvor_cen(i,kt+1)-pvor_cen(i,kt))
97    dyn_tropopause(i) = pplay(i,kt+1)*(pplay(i,kt)/pplay(i,kt+1))**al
98    !--- THETA=380K IN THE TROPICAL REGION
99    al = (tpot_cen(i,kp+1)-th0)/(tpot_cen(i,kp+1)-tpot_cen(i,kp))
100    dyn_tropopause(i) = MAX( pplay(i,kp+1)*(pplay(i,kp)/pplay(i,kp+1))**al,    &
101                            dyn_tropopause(i) )
102  END DO
103
104END FUNCTION dyn_tropopause
105
106
107!-------------------------------------------------------------------------------
108!
109FUNCTION smooth(v,w)
110!
111!-------------------------------------------------------------------------------
112! Arguments:
113  REAL, INTENT(IN)         :: v(:), w(:)
114  REAL, DIMENSION(SIZE(v)) :: smooth
115!-------------------------------------------------------------------------------
116! Local variables:
117  INTEGER :: nv, nw, k, kb, ke, lb, le
118!-------------------------------------------------------------------------------
119  nv=SIZE(v); nw=(SIZE(w)-1)/2
120  DO k=1,nv
121    kb=MAX(k-nw,1 ); lb=MAX(2+nw   -k,1)
122    ke=MIN(k+nw,nv); le=MIN(1+nw+nv-k,1+2*nw)
123    smooth(k)=SUM(v(kb:ke)*w(lb:le))/SUM(w(lb:le))
124  END DO
125
126END FUNCTION smooth
127!
128!-------------------------------------------------------------------------------
129
130END MODULE tropopause_m
Note: See TracBrowser for help on using the repository browser.