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

Last change on this file since 2971 was 2971, checked in by dcugnet, 7 years ago

Control outputs for debug are removed.

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