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