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

Last change on this file since 2930 was 2788, checked in by dcugnet, 8 years ago

Changes in ce0l about the way ozone forcing files are generated:

1) 3D raw input files "climoz.nc" are now handled.
2) Default behaviour is now to let the gcm interpolate in time online.

This helps to avoid huge forcing files (in particular for 3D fields).
In this case, the output files "climoz_LMDZ.nc" all have 14 records:

  • records 2-13 are obtained with records 1-12 of "climoz.nc".
  • records 1 and 14 are obtained respectively with:
    • record 12 of "climoz_m.nc" if available, of "climoz.nc" otherwise.
    • record 1 of "climoz_p.nc" if available, of "climoz.nc" otherwise.

3) If ok_daily_climoz key is TRUE, the time interpolation (one record

a day) is forced, using the 14 records described below.
This now depends on the calendar (it was on a 360 days basis only).

Changes in the gcm about the way zone forcing files are read/interpolated:

1) 3D horizontally interpolated "climoz_LMDZ.nc" files are now handled.
2) Daily files (already interpolated in time) are still handled, but their

number of records must match the expected number of days, that depends
on the calendar (records step is no longer 1/360 year).

3) 14 records monthly files are now handled (and prefered). This reduces

the I/O to a minimum and the aditional computational cost is low (simple
online linear time interpolation).

4) If adjust_tropopause key is TRUE, the input fields are stretched using

following method:

  • LMDZ dynamical tropopause is detected: Ptrop_lmdz = MAX ( P(Potential Vorticity==2PVU), P(theta==380K) )
  • file chemical tropopause is detected: Ptrop_file = P( tro3 == o3t ), where:

o3t = 91. + 28. * SIN(PI*(month-2)/6) (ppbV)

This formula comes from Thouret & al., ACP 6, 1033-1051, 2006.
The second term of the expression is multiplied by TANH(lat_deg/20.)
to account for latitude dependency.

  • File profile is streched in a +/- 5kms zone around the mean tropopause to ensure resulting tropopause matches the one of LMDZ. See procedure regr_pr_time_av for more details.
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.