source: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/tropopause_m.F90 @ 3698

Last change on this file since 3698 was 3141, checked in by dcugnet, 7 years ago
  • Fix: Sig_bot was used before being initialized in "regr_pr_time_av".
  • Add an optional output argument to the "dyn_tropopause" function: the table of the last tropospheric layer in each column.
  • A filter for irrealistic tropopause pressure (lower than 80hPa or greater than 400hPa) is included in "dyn_tropopause". This is safer if this function is not used together with "regr_pr_time_av" (which contains already this filter), which is the case for the LMDZ version with coupled REPROBUS and INCA chemistry.
File size: 5.9 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, itrop, 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  INTEGER, INTENT(OUT), OPTIONAL :: itrop(klon) !--- Last tropospheric layer idx
29  REAL,    INTENT(IN),  OPTIONAL :: thet0, pvor0
30!-------------------------------------------------------------------------------
31! Local variables:
32  include "YOMCST.h"
33  REAL, PARAMETER :: DynPTrMin =8.E+3 !--- Thresholds for minimum and maximum
34  REAL, PARAMETER :: DynPTrMax =4.E+4 !    dynamical tropopause pressure (Pa).
35  CHARACTER(LEN=80)  :: sub
36  INTEGER :: i, k, kb, kt, kp, ib, ie, nw
37  REAL    :: al, th0, pv0
38  REAL,    DIMENSION(klon,klev) :: tpot_cen, tpot_edg, pvor_cen
39  REAL,    PARAMETER :: sg0=0.75  !--- Start level for PV=cte search loop
40  INTEGER, PARAMETER :: nadj=3    !--- Adjacent levs nb for thresholds detection
41  REAL,    PARAMETER :: w(5)=[0.1,0.25,0.3,0.25,0.1] !--- Vertical smoothing
42  INTEGER, SAVE :: k0
43  LOGICAL, SAVE :: first=.TRUE.
44!$OMP THREADPRIVATE(k0,first)
45!-------------------------------------------------------------------------------
46  sub='dyn_tropopause'
47  CALL assert(SIZE(t ,1)==klon, TRIM(sub)//" t klon")
48  CALL assert(SIZE(t ,2)==klev, TRIM(sub)//" t klev")
49  CALL assert(SIZE(ts,1)==klon, TRIM(sub)//" ts klon")
50  CALL assert(SHAPE(paprs)==[klon,klev+1],TRIM(sub)//" paprs shape")
51  CALL assert(SHAPE(pplay)==[klon,klev  ],TRIM(sub)//" pplay shape")
52  CALL assert(SHAPE(rot)  ==[klon,klev  ],TRIM(sub)//" rot shape")
53
54  !--- DEFAULT THRESHOLDS
55  th0=380.; IF(PRESENT(thet0)) th0=thet0   !--- In kelvins
56  pv0=  2.; IF(PRESENT(pvor0)) pv0=pvor0   !--- In PVU
57  IF(first) THEN
58    DO k0=1,klev; IF(aps(k0)/preff+bps(k0)<sg0) EXIT; END DO; first=.FALSE.
59  END IF
60
61  !--- POTENTIAL TEMPERATURE AT CELLS CENTERS AND INTERFACES
62  DO i = 1,klon
63    tpot_cen(i,1) = t(i,1)*(preff/pplay(i,1))**RKAPPA
64    tpot_edg(i,1) = ts(i) *(preff/paprs(i,1))**RKAPPA
65    DO k=2,klev
66      al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
67      tpot_cen(i,k) =  t(i,k)                        *(preff/pplay(i,k))**RKAPPA
68      tpot_edg(i,k) = (t(i,k-1)+al*(t(i,k)-t(i,k-1)))*(preff/paprs(i,k))**RKAPPA
69      !--- FORCE QUANTITIES TO BE GROWING
70      IF(tpot_edg(i,k)<tpot_edg(i,k-1)) tpot_edg(i,k)=tpot_edg(i,k-1)+1.E-5
71      IF(tpot_cen(i,k)<tpot_cen(i,k-1)) tpot_cen(i,k)=tpot_cen(i,k-1)+1.E-5
72    END DO
73    !--- VERTICAL SMOOTHING
74    tpot_cen(i,:)=smooth(tpot_cen(i,:),w)
75    tpot_edg(i,:)=smooth(tpot_edg(i,:),w)
76  END DO
77
78  !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer)
79  DO i = 1, klon
80    DO k= 1, klev-1
81      pvor_cen(i,k)=-1.E6*RG*(rot(i,k)+2.*ROMEGA*SIN(latitude_deg(i)*RPI/180.))&
82                   * (tpot_edg(i,k+1)-tpot_edg(i,k)) / (paprs(i,k+1)-paprs(i,k))
83    END DO
84    !--- VERTICAL SMOOTHING
85    pvor_cen(i,1:klev-1)=smooth(pvor_cen(i,1:klev-1),w)
86  END DO
87
88  !--- LOCATE TROPOPAUSE: LOWEST POINT BETWEEN THETA=380K AND PV=2PVU SURFACES.
89  DO i = 1, klon
90    !--- UPPER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM TOP
91    DO kt=klev-1,1,-1; IF(ALL(ABS(pvor_cen(i,kt-nadj:kt))<=pv0)) EXIT; END DO
92    !--- LOWER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM BOTTOM
93    DO kb=k0,klev-1;   IF(ALL(ABS(pvor_cen(i,kb:kb+nadj))> pv0)) EXIT; END DO; kb=kb-1
94    !--- ISO-THETA POINT: THETA=380K       STARTING FROM TOP
95    DO kp=klev-1,1,-1; IF(ALL(ABS(tpot_cen(i,kp-nadj:kp))<=th0)) EXIT; END DO
96    !--- CHOOSE BETWEEN LOWER AND UPPER TROPOPAUSE
97    IF(2*COUNT(ABS(pvor_cen(i,kb:kt))>pv0)>kt-kb+1) kt=kb
98    !--- PV-DEFINED TROPOPAUSE
99    al = (ABS(pvor_cen(i,kt+1))-pv0)/ABS(pvor_cen(i,kt+1)-pvor_cen(i,kt))
100    dyn_tropopause(i) = pplay(i,kt+1)*(pplay(i,kt)/pplay(i,kt+1))**al
101    !--- THETA=380K IN THE TROPICAL REGION
102    al = (tpot_cen(i,kp+1)-th0)/(tpot_cen(i,kp+1)-tpot_cen(i,kp))
103    dyn_tropopause(i) = MAX( pplay(i,kp+1)*(pplay(i,kp)/pplay(i,kp+1))**al,    &
104                            dyn_tropopause(i) )
105    !--- UNREALISTIC VALUES DETECTION
106    IF(dyn_tropopause(i)<DynPTrMin.OR.dyn_tropopause(i)>DynPTrMax) THEN
107      dyn_tropopause(i)=MIN(MAX(dyn_tropopause(i),DynPTrMax),DynPTrMin)
108      DO kt=1,klev-1; IF(pplay(i,kt+1)>dyn_tropopause(i)) EXIT; END DO; kp=kt
109    END IF
110    !--- LAST TROPOSPHERIC LAYER INDEX NEEDED
111    IF(PRESENT(itrop)) itrop(i)=MAX(kt,kp)
112  END DO
113
114END FUNCTION dyn_tropopause
115
116
117!-------------------------------------------------------------------------------
118!
119FUNCTION smooth(v,w)
120!
121!-------------------------------------------------------------------------------
122! Arguments:
123  REAL, INTENT(IN)         :: v(:), w(:)
124  REAL, DIMENSION(SIZE(v)) :: smooth
125!-------------------------------------------------------------------------------
126! Local variables:
127  INTEGER :: nv, nw, k, kb, ke, lb, le
128!-------------------------------------------------------------------------------
129  nv=SIZE(v); nw=(SIZE(w)-1)/2
130  DO k=1,nv
131    kb=MAX(k-nw,1 ); lb=MAX(2+nw   -k,1)
132    ke=MIN(k+nw,nv); le=MIN(1+nw+nv-k,1+2*nw)
133    smooth(k)=SUM(v(kb:ke)*w(lb:le))/SUM(w(lb:le))
134  END DO
135
136END FUNCTION smooth
137!
138!-------------------------------------------------------------------------------
139
140END MODULE tropopause_m
Note: See TracBrowser for help on using the repository browser.