source: LMDZ6/branches/LMDZ-COSP/libf/phylmd/tropopause_m.F90 @ 5409

Last change on this file since 5409 was 3666, checked in by lfalletti, 5 years ago

Adding changes for Reprobus

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