| 1 | subroutine gwprofil |
|---|
| 2 | * ( nlon, nlev |
|---|
| 3 | * , kgwd ,kdx , ktest |
|---|
| 4 | * , kkcrit, kkcrith, kcrit , kkenvh, kknu,kknu2 |
|---|
| 5 | * , kkbreak |
|---|
| 6 | * , paphm1, prho , pstab , ptfr , pvph , pri , ptau |
|---|
| 7 | * , pdmod , pnu , psig ,pgamma, pstd, ppic,pval) |
|---|
| 8 | |
|---|
| 9 | C**** *gwprofil* |
|---|
| 10 | C |
|---|
| 11 | C purpose. |
|---|
| 12 | C -------- |
|---|
| 13 | C |
|---|
| 14 | C** interface. |
|---|
| 15 | C ---------- |
|---|
| 16 | C from *gwdrag* |
|---|
| 17 | C |
|---|
| 18 | C explicit arguments : |
|---|
| 19 | C -------------------- |
|---|
| 20 | C ==== inputs === |
|---|
| 21 | C |
|---|
| 22 | C ==== outputs === |
|---|
| 23 | C |
|---|
| 24 | C implicit arguments : none |
|---|
| 25 | C -------------------- |
|---|
| 26 | C |
|---|
| 27 | C method: |
|---|
| 28 | C ------- |
|---|
| 29 | C the stress profile for gravity waves is computed as follows: |
|---|
| 30 | C it decreases linearly with heights from the ground |
|---|
| 31 | C to the low-level indicated by kkcrith, |
|---|
| 32 | C to simulates lee waves or |
|---|
| 33 | C low-level gravity wave breaking. |
|---|
| 34 | C above it is constant, except when the waves encounter a critical |
|---|
| 35 | C level (kcrit) or when they break. |
|---|
| 36 | C The stress is also uniformly distributed above the level |
|---|
| 37 | C ntop. |
|---|
| 38 | C |
|---|
| 39 | use dimphy |
|---|
| 40 | use YOMCST_mod |
|---|
| 41 | IMPLICIT NONE |
|---|
| 42 | |
|---|
| 43 | C#include "YOMCST.h" |
|---|
| 44 | #include "YOEGWD.h" |
|---|
| 45 | |
|---|
| 46 | C----------------------------------------------------------------------- |
|---|
| 47 | C |
|---|
| 48 | C* 0.1 ARGUMENTS |
|---|
| 49 | C --------- |
|---|
| 50 | C |
|---|
| 51 | integer nlon,nlev,kgwd |
|---|
| 52 | integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon) |
|---|
| 53 | * ,kdx(nlon),ktest(nlon) |
|---|
| 54 | * ,kkenvh(nlon),kknu(nlon),kknu2(nlon),kkbreak(nlon) |
|---|
| 55 | C |
|---|
| 56 | real paphm1(nlon,nlev+1), pstab(nlon,nlev+1), |
|---|
| 57 | * prho (nlon,nlev+1), pvph (nlon,nlev+1), |
|---|
| 58 | * pri (nlon,nlev+1), ptfr (nlon), ptau(nlon,nlev+1) |
|---|
| 59 | |
|---|
| 60 | real pdmod (nlon) , pnu (nlon) , psig(nlon), |
|---|
| 61 | * pgamma(nlon) , pstd(nlon) , ppic(nlon), pval(nlon) |
|---|
| 62 | |
|---|
| 63 | C----------------------------------------------------------------------- |
|---|
| 64 | C |
|---|
| 65 | C* 0.2 local arrays |
|---|
| 66 | C ------------ |
|---|
| 67 | C |
|---|
| 68 | integer jl,jk |
|---|
| 69 | real zsqr,zalfa,zriw,zdel,zb,zalpha,zdz2n,zdelp,zdelpt |
|---|
| 70 | |
|---|
| 71 | real zdz2 (klon,klev) , znorm(klon) , zoro(klon) |
|---|
| 72 | real ztau (klon,klev+1) |
|---|
| 73 | C |
|---|
| 74 | C----------------------------------------------------------------------- |
|---|
| 75 | C |
|---|
| 76 | C* 1. INITIALIZATION |
|---|
| 77 | C -------------- |
|---|
| 78 | C |
|---|
| 79 | C print *,' entree gwprofil' |
|---|
| 80 | 100 CONTINUE |
|---|
| 81 | C |
|---|
| 82 | C |
|---|
| 83 | C* COMPUTATIONAL CONSTANTS. |
|---|
| 84 | C ------------- ---------- |
|---|
| 85 | C |
|---|
| 86 | do 400 jl=kidia,kfdia |
|---|
| 87 | if(ktest(jl).eq.1)then |
|---|
| 88 | zoro(jl)=psig(jl)*pdmod(jl)/4./pstd(jl) |
|---|
| 89 | ztau(jl,klev+1)=ptau(jl,klev+1) |
|---|
| 90 | c print *,jl,ptau(jl,klev+1) |
|---|
| 91 | ztau(jl,kkcrith(jl))=grahilo*ptau(jl,klev+1) |
|---|
| 92 | endif |
|---|
| 93 | 400 continue |
|---|
| 94 | |
|---|
| 95 | C |
|---|
| 96 | do 430 jk=klev+1,1,-1 |
|---|
| 97 | C |
|---|
| 98 | C |
|---|
| 99 | C* 4.1 constant shear stress until top of the |
|---|
| 100 | C low-level breaking/trapped layer |
|---|
| 101 | 410 CONTINUE |
|---|
| 102 | C |
|---|
| 103 | do 411 jl=kidia,kfdia |
|---|
| 104 | if(ktest(jl).eq.1)then |
|---|
| 105 | if(jk.gt.kkcrith(jl)) then |
|---|
| 106 | zdelp=paphm1(jl,jk)-paphm1(jl,klev+1) |
|---|
| 107 | zdelpt=paphm1(jl,kkcrith(jl))-paphm1(jl,klev+1) |
|---|
| 108 | ptau(jl,jk)=ztau(jl,klev+1)+zdelp/zdelpt* |
|---|
| 109 | c (ztau(jl,kkcrith(jl))-ztau(jl,klev+1)) |
|---|
| 110 | else |
|---|
| 111 | ptau(jl,jk)=ztau(jl,kkcrith(jl)) |
|---|
| 112 | endif |
|---|
| 113 | endif |
|---|
| 114 | 411 continue |
|---|
| 115 | C |
|---|
| 116 | C* 4.15 constant shear stress until the top of the |
|---|
| 117 | C low level flow layer. |
|---|
| 118 | 415 continue |
|---|
| 119 | C |
|---|
| 120 | C |
|---|
| 121 | C* 4.2 wave displacement at next level. |
|---|
| 122 | C |
|---|
| 123 | 420 continue |
|---|
| 124 | C |
|---|
| 125 | 430 continue |
|---|
| 126 | |
|---|
| 127 | C |
|---|
| 128 | C* 4.4 wave richardson number, new wave displacement |
|---|
| 129 | C* and stress: breaking evaluation and critical |
|---|
| 130 | C level |
|---|
| 131 | C |
|---|
| 132 | |
|---|
| 133 | |
|---|
| 134 | do 440 jk=klev,1,-1 |
|---|
| 135 | |
|---|
| 136 | do 441 jl=kidia,kfdia |
|---|
| 137 | if(ktest(jl).eq.1)then |
|---|
| 138 | znorm(jl)=prho(jl,jk)*sqrt(pstab(jl,jk))*pvph(jl,jk) |
|---|
| 139 | zdz2(jl,jk)=ptau(jl,jk)/amax1(znorm(jl),gssec)/zoro(jl) |
|---|
| 140 | endif |
|---|
| 141 | 441 continue |
|---|
| 142 | |
|---|
| 143 | do 442 jl=kidia,kfdia |
|---|
| 144 | if(ktest(jl).eq.1)then |
|---|
| 145 | if(jk.lt.kkcrith(jl)) then |
|---|
| 146 | if((ptau(jl,jk+1).lt.gtsec).or.(jk.le.kcrit(jl))) then |
|---|
| 147 | ptau(jl,jk)=0.0 |
|---|
| 148 | else |
|---|
| 149 | zsqr=sqrt(pri(jl,jk)) |
|---|
| 150 | zalfa=sqrt(pstab(jl,jk)*zdz2(jl,jk))/pvph(jl,jk) |
|---|
| 151 | zriw=pri(jl,jk)*(1.-zalfa)/(1+zalfa*zsqr)**2 |
|---|
| 152 | if(zriw.lt.grcrit) then |
|---|
| 153 | c print *,' breaking!!!',ptau(jl,jk),zsqr |
|---|
| 154 | zdel=4./zsqr/grcrit+1./grcrit**2+4./grcrit |
|---|
| 155 | zb=1./grcrit+2./zsqr |
|---|
| 156 | zalpha=0.5*(-zb+sqrt(zdel)) |
|---|
| 157 | zdz2n=(pvph(jl,jk)*zalpha)**2/pstab(jl,jk) |
|---|
| 158 | ptau(jl,jk)=znorm(jl)*zdz2n*zoro(jl) |
|---|
| 159 | endif |
|---|
| 160 | |
|---|
| 161 | ptau(jl,jk)=amin1(ptau(jl,jk),ptau(jl,jk+1)) |
|---|
| 162 | |
|---|
| 163 | endif |
|---|
| 164 | endif |
|---|
| 165 | endif |
|---|
| 166 | 442 continue |
|---|
| 167 | 440 continue |
|---|
| 168 | |
|---|
| 169 | C REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL |
|---|
| 170 | |
|---|
| 171 | do 530 jl=kidia,kfdia |
|---|
| 172 | if(ktest(jl).eq.1)then |
|---|
| 173 | ztau(jl,kkcrith(jl)-1)=ptau(jl,kkcrith(jl)-1) |
|---|
| 174 | ztau(jl,ntop)=ptau(jl,ntop) |
|---|
| 175 | endif |
|---|
| 176 | 530 continue |
|---|
| 177 | |
|---|
| 178 | do 531 jk=1,klev |
|---|
| 179 | |
|---|
| 180 | do 532 jl=kidia,kfdia |
|---|
| 181 | if(ktest(jl).eq.1)then |
|---|
| 182 | |
|---|
| 183 | if(jk.gt.kkcrith(jl)-1)then |
|---|
| 184 | |
|---|
| 185 | zdelp=paphm1(jl,jk)-paphm1(jl,klev+1 ) |
|---|
| 186 | zdelpt=paphm1(jl,kkcrith(jl)-1)-paphm1(jl,klev+1 ) |
|---|
| 187 | ptau(jl,jk)=ztau(jl,klev+1 ) + |
|---|
| 188 | . (ztau(jl,kkcrith(jl)-1)-ztau(jl,klev+1 ) )* |
|---|
| 189 | . zdelp/zdelpt |
|---|
| 190 | |
|---|
| 191 | endif |
|---|
| 192 | endif |
|---|
| 193 | |
|---|
| 194 | 532 continue |
|---|
| 195 | |
|---|
| 196 | C REORGANISATION AT THE MODEL TOP.... |
|---|
| 197 | |
|---|
| 198 | do 533 jl=kidia,kfdia |
|---|
| 199 | if(ktest(jl).eq.1)then |
|---|
| 200 | |
|---|
| 201 | if(jk.lt.ntop)then |
|---|
| 202 | |
|---|
| 203 | zdelp =paphm1(jl,ntop) |
|---|
| 204 | zdelpt=paphm1(jl,jk) |
|---|
| 205 | ptau(jl,jk)=ztau(jl,ntop)*zdelpt/zdelp |
|---|
| 206 | c ptau(jl,jk)=ztau(jl,ntop) |
|---|
| 207 | |
|---|
| 208 | endif |
|---|
| 209 | |
|---|
| 210 | endif |
|---|
| 211 | |
|---|
| 212 | 533 continue |
|---|
| 213 | |
|---|
| 214 | |
|---|
| 215 | 531 continue |
|---|
| 216 | |
|---|
| 217 | c Yo, this is Venus. |
|---|
| 218 | do jl=kidia,kfdia |
|---|
| 219 | do jk=klev,1,-1 |
|---|
| 220 | if(ktest(jl).eq.1)then |
|---|
| 221 | if(jk.lt.kkbreak(jl)) ptau(jl,jk)=0.0 |
|---|
| 222 | endif |
|---|
| 223 | enddo |
|---|
| 224 | enddo |
|---|
| 225 | |
|---|
| 226 | |
|---|
| 227 | ! Venus: resolve waves |
|---|
| 228 | do jk=klev,1,-1 |
|---|
| 229 | do jl=kidia,kfdia |
|---|
| 230 | if(ktest(jl).eq.1)then |
|---|
| 231 | ! if surface stress greater than threshold |
|---|
| 232 | if (ztau(jl,klev+1) .ge. taubs) then |
|---|
| 233 | ! then enforce same stress in the atmosphere up to the predefined level |
|---|
| 234 | if ((jk.gt.levbs)) then |
|---|
| 235 | ptau(jl,jk) = ztau(jl,klev+1) |
|---|
| 236 | ! and zero above |
|---|
| 237 | elseif (jk.le.levbs) then |
|---|
| 238 | ptau(jl,jk) = 0. |
|---|
| 239 | endif |
|---|
| 240 | ! else |
|---|
| 241 | !if (jk.le.klev-1) ptau(jl,jk) = 0. |
|---|
| 242 | ! ptau(jl,jk) = 0. |
|---|
| 243 | endif |
|---|
| 244 | endif |
|---|
| 245 | enddo |
|---|
| 246 | enddo |
|---|
| 247 | |
|---|
| 248 | |
|---|
| 249 | |
|---|
| 250 | 123 format(i4,1x,20(f6.3,1x)) |
|---|
| 251 | |
|---|
| 252 | |
|---|
| 253 | return |
|---|
| 254 | end |
|---|
| 255 | |
|---|
