source: LMDZ5/branches/testing/libf/phylmd/tropopause_m.F90 @ 5003

Last change on this file since 5003 was 2839, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2785:2838 into testing branch

File size: 4.1 KB
Line 
1MODULE tropopause_m
2
3!  USE phys_local_var_mod, ONLY: ptrop => pr_tropopause
4  PRIVATE
5  PUBLIC :: dyn_tropopause
6
7CONTAINS
8
9!-------------------------------------------------------------------------------
10!
11SUBROUTINE dyn_tropopause(t,ts,paprs,pplay,presnivs,rot,ptrop,tpot0,pvor0,pmin0,pmax0)
12!
13!-------------------------------------------------------------------------------
14  USE assert_m,      ONLY: assert
15  USE assert_eq_m,   ONLY: assert_eq
16  USE comvert_mod,   ONLY: preff
17  USE dimphy,        ONLY: klon, klev
18  USE geometry_mod,  ONLY: latitude_deg, longitude_deg
19  USE interpolation, ONLY: locate
20  IMPLICIT NONE
21!-------------------------------------------------------------------------------
22! Arguments:
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)  :: presnivs(:) !--- Cells-centers nominal pressure
28  REAL, INTENT(IN)  ::    rot(:,:) !--- Cells-centers relative vorticity
29  REAL, INTENT(OUT) :: ptrop(klon) !--- Tropopause pressure
30  REAL, INTENT(IN), OPTIONAL :: tpot0, pvor0, pmin0, pmax0
31!-------------------------------------------------------------------------------
32! Local variables:
33  include "YOMCST.h"
34  CHARACTER(LEN=80)  :: sub
35  INTEGER :: it, i, nsrf, k, nh, kmin, kmax
36  REAL    :: p1, p2, tpo0, pvo0, pmn0, pmx0, al
37  REAL, DIMENSION(klon,klev)   :: t_edg, t_pot
38  REAL, DIMENSION(klon,klev-1) :: pvort
39!-------------------------------------------------------------------------------
40  sub='dyn_tropopause'
41  CALL assert(SIZE(t ,1)==klon, TRIM(sub)//" t klon")
42  CALL assert(SIZE(t ,2)==klev, TRIM(sub)//" t klev")
43  CALL assert(SIZE(ts,1)==klon, TRIM(sub)//" ts klon")
44  CALL assert(SIZE(presnivs,1)==klev, TRIM(sub)//" presnivs klev")
45  CALL assert(SHAPE(paprs)==[klon,klev+1],TRIM(sub)//" paprs shape")
46  CALL assert(SHAPE(pplay)==[klon,klev  ],TRIM(sub)//" pplay shape")
47  CALL assert(SHAPE(rot)  ==[klon,klev  ],TRIM(sub)//" rot shape")
48
49  !--- DEFAULT THRESHOLDS
50  tpo0=380.;   IF(PRESENT(tpot0)) tpo0=tpot0  !--- In kelvins
51  pvo0=2.0;    IF(PRESENT(pvor0)) pvo0=pvor0  !--- In PVU
52  pmn0= 8000.; IF(PRESENT(pmin0)) pmn0=pmin0  !--- In pascals
53  pmx0=50000.; IF(PRESENT(pmax0)) pmx0=pmax0  !--- In pascals
54  kmin=klev-locate(presnivs(klev:1:-1),pmx0)+1
55  kmax=klev-locate(presnivs(klev:1:-1),pmn0)+1
56
57  !--- POTENTIAL TEMPERATURE AT CELLS CENTERS
58  DO k= 1, klev
59    DO i = 1, klon
60      t_pot(i,k) = t(i,k)*(preff/pplay(i,k))**RKAPPA
61    END DO
62  END DO
63
64  !--- TEMPERATURE AT CELLS INTERFACES (except in top layer)
65  t_edg(:,1) = ts(:)
66  DO k= 2, klev
67    DO i = 1, klon
68      al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
69      t_edg(i,k) = t(i,k) + al*(t(i,k-1)-t(i,k))
70    END DO
71  END DO
72
73  !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer)
74  DO k= 1, klev-1
75    DO i = 1, klon
76      IF(paprs(i,k)==paprs(i,k+1)) THEN; pvort(i,k)=HUGE(1.); CYCLE; END IF
77      pvort(i,k) = -1.E6 * RG * 2.*ROMEGA*ABS(SIN(latitude_deg(i)*RPI/180.))   &
78                         * ( t_edg(i,k+1)*(preff/paprs(i,k+1) )**RKAPPA        &
79                           - t_edg(i,k  )*(preff/paprs(i,k  ) )**RKAPPA)       &
80                         / ( paprs(i,k+1)  -     paprs(i,k  ) )
81    END DO
82  END DO
83
84  !--- LOCATE TROPOPAUSE
85  ptrop(:)=0.
86  DO i = 1, klon
87    !--- Dynamical tropopause
88    it=kmax; DO WHILE(pvort(i,it)>pvo0.AND.it>=kmin); it=it-1; END DO
89    IF(pvort(i,it+1)/=pvort(i,it).AND.ALL([kmax,kmin-1]/=it) &
90      .AND.ALL([pvort(i,it),pvort(i,it+1)]/=HUGE(1.))) THEN
91      al = (pvort(i,it+1)-pvo0)/(pvort(i,it+1)-pvort(i,it))
92      ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
93    END IF
94    !--- Potential temperature iso-surface
95    it = kmin-1+locate(t_pot(i,kmin:kmax),tpo0)
96    IF(t_pot(i,it+1)/=t_pot(i,it).AND.ALL([kmin-1,kmax]/=it)) THEN
97      al = (t_pot(i,it+1)-tpo0)/(t_pot(i,it+1)-t_pot(i,it))
98      ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
99    END IF
100  END DO
101
102END SUBROUTINE dyn_tropopause
103
104END MODULE tropopause_m
Note: See TracBrowser for help on using the repository browser.