| 1 | MODULE convection |
|---|
| 2 | |
|---|
| 3 | #include "use_logging.h" |
|---|
| 4 | |
|---|
| 5 | IMPLICIT NONE |
|---|
| 6 | |
|---|
| 7 | CONTAINS |
|---|
| 8 | |
|---|
| 9 | PURE SUBROUTINE sigma_levels(ngrid, nlay, i, pplev, ppopsk, sig, dsig, sdsig) |
|---|
| 10 | INTEGER, INTENT(IN) :: ngrid, nlay, i |
|---|
| 11 | REAL, INTENT(IN) :: pplev(ngrid,nlay+1), ppopsk(ngrid,nlay) |
|---|
| 12 | REAL, INTENT(OUT) :: sig(nlay+1), dsig(nlay), sdsig(nlay) |
|---|
| 13 | ! Calcul des niveaux sigma sur cette colonne |
|---|
| 14 | INTEGER :: l |
|---|
| 15 | DO l=1,nlay+1 |
|---|
| 16 | sig(l)=pplev(i,l)/pplev(i,1) |
|---|
| 17 | ENDDO |
|---|
| 18 | DO l=1,nlay |
|---|
| 19 | dsig(l)=sig(l)-sig(l+1) |
|---|
| 20 | sdsig(l)=ppopsk(i,l)*dsig(l) |
|---|
| 21 | ENDDO |
|---|
| 22 | END SUBROUTINE sigma_levels |
|---|
| 23 | |
|---|
| 24 | SUBROUTINE adjust_onecolumn(ngrid, nlay, i, sig, dsig, sdsig, zu2, zv2, zh2) |
|---|
| 25 | INTEGER, INTENT(IN) :: ngrid, nlay, i |
|---|
| 26 | REAL, INTENT(OUT) :: sig(nlay+1), sdsig(nlay), dsig(nlay) |
|---|
| 27 | REAL, INTENT(INOUT) :: zu2(ngrid,nlay), zv2(ngrid,nlay), zh2(ngrid,nlay) |
|---|
| 28 | |
|---|
| 29 | INTEGER :: l,l1,l2 |
|---|
| 30 | LOGICAL :: extend |
|---|
| 31 | REAL :: zhm,zsm,zum,zvm,zalpha |
|---|
| 32 | |
|---|
| 33 | l2=1 |
|---|
| 34 | DO WHILE( l2<nlay ) |
|---|
| 35 | l2 = l2+1 |
|---|
| 36 | ! starting from l1=l2, we shall extend l1..l2 downwards and upwards |
|---|
| 37 | ! until zhm = potential temperature mass-weighted over l1..l2 be |
|---|
| 38 | ! higher than pot. temp. at level l1-1 |
|---|
| 39 | ! and lower than pot. temp. at level l2+1 |
|---|
| 40 | ! for this we incrementally compute zsm = mass of layers l1 .. l2 and zhm |
|---|
| 41 | zsm = sdsig(l2) |
|---|
| 42 | zhm = zh2(i, l2) |
|---|
| 43 | l1 = l2 |
|---|
| 44 | extend = .TRUE. |
|---|
| 45 | DO WHILE(extend) |
|---|
| 46 | ! extend downwards if l1>1 and zhm is lower than layer l1-1 |
|---|
| 47 | extend = .FALSE. |
|---|
| 48 | IF (l1 > 1) THEN |
|---|
| 49 | IF (zhm < zh2(i, l1-1)) extend = .TRUE. |
|---|
| 50 | END IF |
|---|
| 51 | IF(extend) THEN |
|---|
| 52 | l1 = l1-1 |
|---|
| 53 | zsm = zsm + sdsig(l1) |
|---|
| 54 | zhm = zhm + sdsig(l1) * (zh2(i, l1) - zhm) / zsm |
|---|
| 55 | CYCLE |
|---|
| 56 | END IF |
|---|
| 57 | ! extend upwards if l2<nlay and zhm is higher than layer l2+1 |
|---|
| 58 | extend=.FALSE. |
|---|
| 59 | IF (l2 < nlay) THEN |
|---|
| 60 | IF (zh2(i, l2+1) < zhm) extend=.TRUE. |
|---|
| 61 | END IF |
|---|
| 62 | IF(extend) THEN |
|---|
| 63 | l2 = l2+1 |
|---|
| 64 | zsm = zsm + sdsig(l2) |
|---|
| 65 | zhm = zhm + sdsig(l2) * (zh2(i, l2) - zhm) / zsm |
|---|
| 66 | END IF |
|---|
| 67 | END DO |
|---|
| 68 | |
|---|
| 69 | ! at this point the profile l1-1 (l1 ...l2) l2+1 is stable |
|---|
| 70 | |
|---|
| 71 | IF(l1<l2) THEN |
|---|
| 72 | WRITELOG(*,*) 'Mixing between layers ',l1,l2 |
|---|
| 73 | LOG_DBG('convadj') |
|---|
| 74 | ! perform the mixing of layers l1...l2 : |
|---|
| 75 | ! * potential temperature is fully mixed, ie replaced by mass-weighted average |
|---|
| 76 | ! * momentum u,v is mixed with weight zalpha, i.e. u:=zalpha*mean(u)+(1-zalpha)*u |
|---|
| 77 | ! where zalpha = sum( abs(theta-mean(theta))*dsig) / mean(theta)*sum(dsig) |
|---|
| 78 | ! for large deviations of theta from its mean, this formula could in principe yield zalpha>1. |
|---|
| 79 | zalpha=0. |
|---|
| 80 | zum=0. |
|---|
| 81 | zvm=0. |
|---|
| 82 | DO l = l1, l2 |
|---|
| 83 | zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l) |
|---|
| 84 | zh2(i, l) = zhm |
|---|
| 85 | ! we must use compute zum, zvm from zu2, zv2 to conserve momentum (see below) |
|---|
| 86 | zum=zum+dsig(l)*zu2(i,l) |
|---|
| 87 | zvm=zvm+dsig(l)*zv2(i,l) |
|---|
| 88 | END DO |
|---|
| 89 | zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1))) |
|---|
| 90 | zum=zum/(sig(l1)-sig(l2+1)) |
|---|
| 91 | zvm=zvm/(sig(l1)-sig(l2+1)) |
|---|
| 92 | IF(zalpha.GT.1.) THEN |
|---|
| 93 | WRITELOG(*,*) 'zalpha=',zalpha |
|---|
| 94 | CALL log_gridpoint(i) |
|---|
| 95 | LOG_WARN('convadj') |
|---|
| 96 | STOP |
|---|
| 97 | END IF |
|---|
| 98 | IF(zalpha.LT.1.e-5) zalpha=1.e-5 |
|---|
| 99 | ! zum, zvm are mass-weighted averages of zu2, zv2 => mixing preserves total momentum |
|---|
| 100 | DO l=l1,l2 |
|---|
| 101 | zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l)) |
|---|
| 102 | zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l)) |
|---|
| 103 | ENDDO |
|---|
| 104 | END IF |
|---|
| 105 | |
|---|
| 106 | END DO |
|---|
| 107 | |
|---|
| 108 | END SUBROUTINE adjust_onecolumn |
|---|
| 109 | |
|---|
| 110 | SUBROUTINE convadj(ngrid,nlay,ptimestep, & |
|---|
| 111 | & pplay,pplev,ppopsk, & |
|---|
| 112 | & pu,pv,ph, & |
|---|
| 113 | & pdufi,pdvfi,pdhfi, & |
|---|
| 114 | & pduadj,pdvadj,pdhadj) |
|---|
| 115 | |
|---|
| 116 | !======================================================================= |
|---|
| 117 | ! |
|---|
| 118 | ! dry convective adjustment |
|---|
| 119 | ! h = pot. temp. , static instability if ph(above)<ph(below) |
|---|
| 120 | ! tendencies pdfX are first added to profiles pX, X=u,v,h, yieldig pX2 |
|---|
| 121 | ! adjustment is performed on profiles pX2 |
|---|
| 122 | ! then the difference between adjusted and non-adujsted profiles |
|---|
| 123 | ! is converted back as tendencies |
|---|
| 124 | ! |
|---|
| 125 | !======================================================================= |
|---|
| 126 | |
|---|
| 127 | INTEGER, INTENT(IN) :: ngrid,nlay |
|---|
| 128 | REAL, INTENT(IN) :: ptimestep |
|---|
| 129 | REAL, INTENT(IN) :: pplay(ngrid,nlay), pplev(ngrid,nlay+1), ppopsk(ngrid,nlay) |
|---|
| 130 | REAL, INTENT(IN) :: ph(ngrid,nlay), pu(ngrid,nlay), pv(ngrid,nlay), & |
|---|
| 131 | & pdhfi(ngrid,nlay), pdufi(ngrid,nlay), pdvfi(ngrid,nlay) |
|---|
| 132 | REAL, INTENT(OUT) :: pdhadj(ngrid,nlay), pduadj(ngrid,nlay), pdvadj(ngrid,nlay) |
|---|
| 133 | |
|---|
| 134 | ! local: |
|---|
| 135 | REAL sig(nlay+1),sdsig(nlay),dsig(nlay) |
|---|
| 136 | REAL zu(ngrid,nlay), zv(ngrid,nlay), zh(ngrid,nlay) |
|---|
| 137 | REAL zu2(ngrid,nlay), zv2(ngrid,nlay), zh2(ngrid,nlay) |
|---|
| 138 | LOGICAL vtest(ngrid) |
|---|
| 139 | INTEGER jadrs(ngrid) |
|---|
| 140 | INTEGER jcnt, ig, l, jj |
|---|
| 141 | |
|---|
| 142 | ! Add physics tendencies to u,v,h |
|---|
| 143 | DO l=1,nlay |
|---|
| 144 | DO ig=1,ngrid |
|---|
| 145 | zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep |
|---|
| 146 | zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep |
|---|
| 147 | zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep |
|---|
| 148 | END DO |
|---|
| 149 | END DO |
|---|
| 150 | zu2(:,:)=zu(:,:) |
|---|
| 151 | zv2(:,:)=zv(:,:) |
|---|
| 152 | zh2(:,:)=zh(:,:) |
|---|
| 153 | |
|---|
| 154 | !----------------------------------------------------------------------- |
|---|
| 155 | ! detection des profils a modifier: |
|---|
| 156 | ! --------------------------------- |
|---|
| 157 | ! si le profil est a modifier |
|---|
| 158 | ! (i.e. ph(niv_sup) < ph(niv_inf) ) |
|---|
| 159 | ! alors le tableau "vtest" est mis a .TRUE. ; |
|---|
| 160 | ! sinon, il reste a sa valeur initiale (.FALSE.) |
|---|
| 161 | ! cette operation est vectorisable |
|---|
| 162 | |
|---|
| 163 | vtest(:)=.FALSE. |
|---|
| 164 | DO l=2,nlay |
|---|
| 165 | DO ig=1,ngrid |
|---|
| 166 | IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE. |
|---|
| 167 | END DO |
|---|
| 168 | END DO |
|---|
| 169 | |
|---|
| 170 | !---------------------------------- |
|---|
| 171 | ! Ajustement des profils instables |
|---|
| 172 | ! -------------------------------- |
|---|
| 173 | |
|---|
| 174 | DO ig=1,ngrid |
|---|
| 175 | IF(vtest(ig)) THEN |
|---|
| 176 | CALL sigma_levels(ngrid, nlay, ig, pplev, ppopsk, sig, dsig, sdsig) |
|---|
| 177 | CALL adjust_onecolumn(ngrid, nlay, ig, sig, dsig, sdsig, zu2, zv2, zh2) |
|---|
| 178 | ENDIF |
|---|
| 179 | END DO |
|---|
| 180 | |
|---|
| 181 | DO l=1,nlay |
|---|
| 182 | DO ig=1,ngrid |
|---|
| 183 | pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep |
|---|
| 184 | pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep |
|---|
| 185 | pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep |
|---|
| 186 | END DO |
|---|
| 187 | END DO |
|---|
| 188 | |
|---|
| 189 | END SUBROUTINE convadj |
|---|
| 190 | |
|---|
| 191 | END MODULE convection |
|---|
