| [814] | 1 | ! |
|---|
| 2 | ! $Header$ |
|---|
| 3 | ! |
|---|
| 4 | MODULE climb_wind_mod |
|---|
| 5 | ! |
|---|
| 6 | ! Module to solve the verctical diffusion of the wind components "u" and "v". |
|---|
| 7 | ! |
|---|
| 8 | USE dimphy |
|---|
| 9 | |
|---|
| 10 | IMPLICIT NONE |
|---|
| 11 | |
|---|
| 12 | SAVE |
|---|
| 13 | PRIVATE |
|---|
| 14 | |
|---|
| 15 | REAL, DIMENSION(:), ALLOCATABLE :: alf1, alf2 |
|---|
| 16 | !$OMP THREADPRIVATE(alf1,alf2) |
|---|
| 17 | REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefm |
|---|
| 18 | !$OMP THREADPRIVATE(Kcoefm) |
|---|
| 19 | REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_U, Dcoef_U |
|---|
| 20 | !$OMP THREADPRIVATE(Ccoef_U, Dcoef_U) |
|---|
| 21 | REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_V, Dcoef_V |
|---|
| 22 | !$OMP THREADPRIVATE(Ccoef_V, Dcoef_V) |
|---|
| 23 | LOGICAL :: firstcall=.TRUE. |
|---|
| 24 | !$OMP THREADPRIVATE(firstcall) |
|---|
| 25 | |
|---|
| 26 | |
|---|
| 27 | PUBLIC :: climb_wind_down, calcul_wind_flux, climb_wind_up |
|---|
| 28 | |
|---|
| 29 | CONTAINS |
|---|
| 30 | ! |
|---|
| 31 | !**************************************************************************************** |
|---|
| 32 | ! |
|---|
| 33 | SUBROUTINE climb_wind_init |
|---|
| 34 | |
|---|
| 35 | INTEGER :: ierr |
|---|
| 36 | CHARACTER(len = 20) :: modname = 'climb_wind_init' |
|---|
| 37 | |
|---|
| 38 | !**************************************************************************************** |
|---|
| 39 | ! Allocation of global module variables |
|---|
| 40 | ! |
|---|
| 41 | !**************************************************************************************** |
|---|
| 42 | |
|---|
| 43 | ALLOCATE(alf1(klon), stat=ierr) |
|---|
| 44 | IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate alf2',1) |
|---|
| 45 | |
|---|
| 46 | ALLOCATE(alf2(klon), stat=ierr) |
|---|
| 47 | IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate alf2',1) |
|---|
| 48 | |
|---|
| 49 | ALLOCATE(Kcoefm(klon,klev), stat=ierr) |
|---|
| 50 | IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate Kcoefm',1) |
|---|
| 51 | |
|---|
| 52 | ALLOCATE(Ccoef_U(klon,klev), stat=ierr) |
|---|
| 53 | IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate Ccoef_U',1) |
|---|
| 54 | |
|---|
| 55 | ALLOCATE(Dcoef_U(klon,klev), stat=ierr) |
|---|
| 56 | IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocation Dcoef_U',1) |
|---|
| 57 | |
|---|
| 58 | ALLOCATE(Ccoef_V(klon,klev), stat=ierr) |
|---|
| 59 | IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocation Ccoef_V',1) |
|---|
| 60 | |
|---|
| 61 | ALLOCATE(Dcoef_V(klon,klev), stat=ierr) |
|---|
| 62 | IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocation Dcoef_V',1) |
|---|
| 63 | |
|---|
| 64 | firstcall=.FALSE. |
|---|
| 65 | |
|---|
| 66 | END SUBROUTINE climb_wind_init |
|---|
| 67 | ! |
|---|
| 68 | !**************************************************************************************** |
|---|
| 69 | ! |
|---|
| 70 | SUBROUTINE climb_wind_down(knon, dtime, coef_in, pplay, paprs, temp, delp, u_old, v_old) |
|---|
| 71 | ! |
|---|
| 72 | ! This routine calculates for the wind components u and v, |
|---|
| 73 | ! recursivly the coefficients C and D in equation |
|---|
| 74 | ! X(k) = C(k) + D(k)*X(k-1), X=[u,v], k=[1,klev] is the vertical layer. |
|---|
| 75 | ! |
|---|
| 76 | ! |
|---|
| 77 | INCLUDE "YOMCST.h" |
|---|
| 78 | ! Input arguments |
|---|
| 79 | !**************************************************************************************** |
|---|
| 80 | INTEGER, INTENT(IN) :: knon |
|---|
| 81 | REAL, INTENT(IN) :: dtime |
|---|
| 82 | REAL, DIMENSION(klon,klev), INTENT(IN) :: coef_in |
|---|
| 83 | REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pres au milieu de couche (Pa) |
|---|
| 84 | REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche (Pa) |
|---|
| 85 | REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature |
|---|
| 86 | REAL, DIMENSION(klon,klev), INTENT(IN) :: delp |
|---|
| 87 | REAL, DIMENSION(klon,klev), INTENT(IN) :: u_old |
|---|
| 88 | REAL, DIMENSION(klon,klev), INTENT(IN) :: v_old |
|---|
| 89 | |
|---|
| 90 | ! Local variables |
|---|
| 91 | !**************************************************************************************** |
|---|
| 92 | REAL, DIMENSION(klon) :: u1lay, v1lay |
|---|
| 93 | INTEGER :: k, i |
|---|
| 94 | |
|---|
| 95 | |
|---|
| 96 | !**************************************************************************************** |
|---|
| 97 | ! Initialize module |
|---|
| 98 | IF (firstcall) CALL climb_wind_init |
|---|
| 99 | |
|---|
| 100 | !**************************************************************************************** |
|---|
| 101 | ! Calculate the coefficients C and D in : u(k) = C(k) + D(k)*u(k-1) |
|---|
| 102 | ! |
|---|
| 103 | !**************************************************************************************** |
|---|
| 104 | |
|---|
| 105 | ! - Define alpha (alf1 and alf2) |
|---|
| 106 | alf1(:) = 1.0 |
|---|
| 107 | alf2(:) = 1.0 - alf1(:) |
|---|
| 108 | |
|---|
| 109 | ! - Calculte the wind components for the first layer |
|---|
| 110 | u1lay(1:knon) = u_old(1:knon,1)*alf1(1:knon) + u_old(1:knon,2)*alf2(1:knon) |
|---|
| 111 | v1lay(1:knon) = v_old(1:knon,1)*alf1(1:knon) + v_old(1:knon,2)*alf2(1:knon) |
|---|
| 112 | |
|---|
| 113 | ! - Calculate K |
|---|
| 114 | Kcoefm(:,:) = 0.0 |
|---|
| 115 | DO i = 1, knon |
|---|
| 116 | Kcoefm(i,1) = coef_in(i,1) * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))* & |
|---|
| 117 | pplay(i,1)/(RD*temp(i,1)) |
|---|
| 118 | Kcoefm(i,1) = Kcoefm(i,1) * dtime*RG |
|---|
| 119 | END DO |
|---|
| 120 | |
|---|
| 121 | DO k = 2, klev |
|---|
| 122 | DO i=1,knon |
|---|
| 123 | Kcoefm(i,k) = coef_in(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) & |
|---|
| 124 | *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2 |
|---|
| 125 | Kcoefm(i,k) = Kcoefm(i,k) * dtime*RG |
|---|
| 126 | END DO |
|---|
| 127 | END DO |
|---|
| 128 | |
|---|
| 129 | ! - Calculate the coefficients C and D, component "u" |
|---|
| 130 | CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), & |
|---|
| 131 | u_old(:,:), alf1(:), alf2(:), & |
|---|
| 132 | Ccoef_U(:,:), Dcoef_U(:,:)) |
|---|
| 133 | |
|---|
| 134 | ! - Calculate the coefficients C and D, component "v" |
|---|
| 135 | CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), & |
|---|
| 136 | v_old(:,:), alf1(:), alf2(:), & |
|---|
| 137 | Ccoef_V(:,:), Dcoef_V(:,:)) |
|---|
| 138 | |
|---|
| 139 | END SUBROUTINE climb_wind_down |
|---|
| 140 | ! |
|---|
| 141 | !**************************************************************************************** |
|---|
| 142 | ! |
|---|
| 143 | SUBROUTINE calc_coef(knon, Kcoef, dels, X, alfa1, alfa2, Ccoef, Dcoef) |
|---|
| 144 | ! |
|---|
| 145 | ! Find the coefficients C and D in fonction of alfa, K and dels |
|---|
| 146 | ! |
|---|
| 147 | ! Input arguments |
|---|
| 148 | !**************************************************************************************** |
|---|
| 149 | INTEGER, INTENT(IN) :: knon |
|---|
| 150 | REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef, dels |
|---|
| 151 | REAL, DIMENSION(klon,klev), INTENT(IN) :: X |
|---|
| 152 | REAL, DIMENSION(klon), INTENT(IN) :: alfa1, alfa2 |
|---|
| 153 | |
|---|
| 154 | ! Output arguments |
|---|
| 155 | !**************************************************************************************** |
|---|
| 156 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef, Dcoef |
|---|
| 157 | |
|---|
| 158 | ! local variables |
|---|
| 159 | !**************************************************************************************** |
|---|
| 160 | INTEGER :: k, i |
|---|
| 161 | REAL :: buf |
|---|
| 162 | |
|---|
| 163 | !**************************************************************************************** |
|---|
| 164 | ! |
|---|
| 165 | ! Niveau au sommet, k=klev |
|---|
| 166 | ! |
|---|
| 167 | Ccoef(:,:) = 0.0 |
|---|
| 168 | Dcoef(:,:) = 0.0 |
|---|
| 169 | |
|---|
| 170 | DO i = 1, knon |
|---|
| 171 | buf = dels(i,klev) + Kcoef(i,klev) |
|---|
| 172 | |
|---|
| 173 | Ccoef(i,klev) = X(i,klev)*dels(i,klev)/buf |
|---|
| 174 | Dcoef(i,klev) = Kcoef(i,klev)/buf |
|---|
| 175 | END DO |
|---|
| 176 | |
|---|
| 177 | ! |
|---|
| 178 | ! Niveau (klev-1) <= k <= 2 |
|---|
| 179 | ! |
|---|
| 180 | DO k=(klev-1),2,-1 |
|---|
| 181 | DO i = 1, knon |
|---|
| 182 | buf = dels(i,k) + Kcoef(i,k) + & |
|---|
| 183 | Kcoef(i,k+1)*(1.-Dcoef(i,k+1)) |
|---|
| 184 | |
|---|
| 185 | Ccoef(i,k) = X(i,k)*dels(i,k) + & |
|---|
| 186 | Kcoef(i,k+1)*Ccoef(i,k+1) |
|---|
| 187 | Ccoef(i,k) = Ccoef(i,k)/buf |
|---|
| 188 | Dcoef(i,k) = Kcoef(i,k)/buf |
|---|
| 189 | END DO |
|---|
| 190 | END DO |
|---|
| 191 | |
|---|
| 192 | ! |
|---|
| 193 | ! Niveau k=1 |
|---|
| 194 | ! |
|---|
| 195 | DO i = 1, knon |
|---|
| 196 | buf = dels(i,1) + & |
|---|
| 197 | (alfa1(i) + alfa2(i)*Dcoef(i,2)) * Kcoef(i,1) + & |
|---|
| 198 | (1.-Dcoef(i,2))*Kcoef(i,2) |
|---|
| 199 | |
|---|
| 200 | Ccoef(i,1) = X(i,1)*dels(i,1) + & |
|---|
| 201 | (Kcoef(i,2)-Kcoef(i,1)*alfa2(i)) * Ccoef(i,2) |
|---|
| 202 | Ccoef(i,1) = Ccoef(i,1)/buf |
|---|
| 203 | Dcoef(i,1) = Kcoef(i,1)/buf |
|---|
| 204 | END DO |
|---|
| 205 | |
|---|
| 206 | END SUBROUTINE calc_coef |
|---|
| 207 | ! |
|---|
| 208 | !**************************************************************************************** |
|---|
| 209 | ! |
|---|
| 210 | SUBROUTINE calcul_wind_flux(knon, dtime, & |
|---|
| 211 | flux_u, flux_v) |
|---|
| 212 | |
|---|
| 213 | INCLUDE "YOMCST.h" |
|---|
| 214 | |
|---|
| 215 | ! Input arguments |
|---|
| 216 | !**************************************************************************************** |
|---|
| 217 | INTEGER, INTENT(IN) :: knon |
|---|
| 218 | REAL, INTENT(IN) :: dtime |
|---|
| 219 | |
|---|
| 220 | ! Output arguments |
|---|
| 221 | !**************************************************************************************** |
|---|
| 222 | REAL, DIMENSION(klon), INTENT(OUT) :: flux_u |
|---|
| 223 | REAL, DIMENSION(klon), INTENT(OUT) :: flux_v |
|---|
| 224 | |
|---|
| 225 | ! Local variables |
|---|
| 226 | !**************************************************************************************** |
|---|
| 227 | INTEGER :: i |
|---|
| 228 | REAL, DIMENSION(klon) :: u0, v0 ! u and v at niveau 0 |
|---|
| 229 | REAL, DIMENSION(klon) :: u1, v1 ! u and v at niveau 1 |
|---|
| 230 | REAL, DIMENSION(klon) :: u2, v2 ! u and v at niveau 2 |
|---|
| 231 | |
|---|
| 232 | |
|---|
| 233 | !**************************************************************************************** |
|---|
| 234 | ! Les vents de surface sont supposes nuls |
|---|
| 235 | ! |
|---|
| 236 | !**************************************************************************************** |
|---|
| 237 | u0(:) = 0.0 |
|---|
| 238 | v0(:) = 0.0 |
|---|
| 239 | |
|---|
| 240 | !**************************************************************************************** |
|---|
| 241 | ! On calcule les vents du couhes 1 et 2 recurviement |
|---|
| 242 | ! |
|---|
| 243 | !**************************************************************************************** |
|---|
| 244 | DO i = 1, knon |
|---|
| 245 | u1(i) = Ccoef_U(i,1) + Dcoef_U(i,1)*u0(i) |
|---|
| 246 | v1(i) = Ccoef_V(i,1) + Dcoef_V(i,1)*v0(i) |
|---|
| 247 | u2(i) = Ccoef_U(i,2) + Dcoef_U(i,2)*u1(i) |
|---|
| 248 | v2(i) = Ccoef_V(i,2) + Dcoef_V(i,2)*v1(i) |
|---|
| 249 | END DO |
|---|
| 250 | |
|---|
| 251 | !**************************************************************************************** |
|---|
| 252 | ! On calcule le flux |
|---|
| 253 | ! |
|---|
| 254 | !**************************************************************************************** |
|---|
| 255 | flux_u(:) = 0.0 |
|---|
| 256 | flux_v(:) = 0.0 |
|---|
| 257 | |
|---|
| 258 | DO i=1,knon |
|---|
| 259 | flux_u(i) = Kcoefm(i,1)/RG/dtime * (u1(i)*alf1(i) + u2(i)*alf2(i) - u0(i)) |
|---|
| 260 | flux_v(i) = Kcoefm(i,1)/RG/dtime * (v1(i)*alf1(i) + v2(i)*alf2(i) - v0(i)) |
|---|
| 261 | END DO |
|---|
| 262 | |
|---|
| 263 | END SUBROUTINE calcul_wind_flux |
|---|
| 264 | ! |
|---|
| 265 | !**************************************************************************************** |
|---|
| 266 | ! |
|---|
| 267 | SUBROUTINE climb_wind_up(knon, dtime, u_old, v_old, & |
|---|
| 268 | flx_u_new, flx_v_new, d_u_new, d_v_new) |
|---|
| 269 | ! |
|---|
| 270 | ! Diffuse the wind components from the surface and up to the top layer. Coefficents |
|---|
| 271 | ! C and D are known from before. The values for U and V at surface are supposed to be |
|---|
| 272 | ! zero (this could be modified). |
|---|
| 273 | ! |
|---|
| 274 | ! u(k) = Cu(k) + Du(k)*u(k-1) |
|---|
| 275 | ! v(k) = Cv(k) + Dv(k)*v(k-1) |
|---|
| 276 | ! [1 <= k <= klev] |
|---|
| 277 | ! |
|---|
| 278 | !**************************************************************************************** |
|---|
| 279 | INCLUDE "YOMCST.h" |
|---|
| 280 | |
|---|
| 281 | ! Input arguments |
|---|
| 282 | !**************************************************************************************** |
|---|
| 283 | INTEGER, INTENT(IN) :: knon |
|---|
| 284 | REAL, INTENT(IN) :: dtime |
|---|
| 285 | REAL, DIMENSION(klon,klev), INTENT(IN) :: u_old |
|---|
| 286 | REAL, DIMENSION(klon,klev), INTENT(IN) :: v_old |
|---|
| 287 | |
|---|
| 288 | ! Output arguments |
|---|
| 289 | !**************************************************************************************** |
|---|
| 290 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: flx_u_new, flx_v_new |
|---|
| 291 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_u_new, d_v_new |
|---|
| 292 | |
|---|
| 293 | ! Local variables |
|---|
| 294 | !**************************************************************************************** |
|---|
| 295 | REAL, DIMENSION(klon,klev) :: u_new, v_new |
|---|
| 296 | REAL, DIMENSION(klon) :: u0, v0 |
|---|
| 297 | INTEGER :: k, i |
|---|
| 298 | |
|---|
| 299 | ! |
|---|
| 300 | !**************************************************************************************** |
|---|
| 301 | |
|---|
| 302 | ! Niveau 0 |
|---|
| 303 | u0(1:knon) = 0.0 |
|---|
| 304 | v0(1:knon) = 0.0 |
|---|
| 305 | |
|---|
| 306 | ! Niveau 1 |
|---|
| 307 | DO i = 1, knon |
|---|
| 308 | u_new(i,1) = Ccoef_U(i,1) + Dcoef_U(i,1) * u0(i) |
|---|
| 309 | v_new(i,1) = Ccoef_V(i,1) + Dcoef_V(i,1) * v0(i) |
|---|
| 310 | END DO |
|---|
| 311 | |
|---|
| 312 | ! Niveau 2 jusqu'au sommet klev |
|---|
| 313 | DO k = 2, klev |
|---|
| 314 | DO i=1, knon |
|---|
| 315 | u_new(i,k) = Ccoef_U(i,k) + Dcoef_U(i,k) * u_new(i,k-1) |
|---|
| 316 | v_new(i,k) = Ccoef_V(i,k) + Dcoef_V(i,k) * v_new(i,k-1) |
|---|
| 317 | END DO |
|---|
| 318 | END DO |
|---|
| 319 | |
|---|
| 320 | !**************************************************************************************** |
|---|
| 321 | ! Calcul flux |
|---|
| 322 | ! |
|---|
| 323 | !== flux_u/v est le flux de moment angulaire (positif vers bas) |
|---|
| 324 | !== dont l'unite est: (kg m/s)/(m**2 s) |
|---|
| 325 | ! |
|---|
| 326 | !**************************************************************************************** |
|---|
| 327 | ! |
|---|
| 328 | flx_u_new(:,:) = 0.0 |
|---|
| 329 | flx_v_new(:,:) = 0.0 |
|---|
| 330 | |
|---|
| 331 | ! Niveau 1 |
|---|
| 332 | DO i = 1, knon |
|---|
| 333 | flx_u_new(i,1) = Kcoefm(i,1)/RG/dtime * & |
|---|
| 334 | (u_new(i,1)*alf1(i)+u_new(i,2)*alf2(i) - u0(i)) |
|---|
| 335 | flx_v_new(i,1) = Kcoefm(i,1)/RG/dtime * & |
|---|
| 336 | (v_new(i,1)*alf1(i)+v_new(i,2)*alf2(i) - v0(i)) |
|---|
| 337 | END DO |
|---|
| 338 | |
|---|
| 339 | ! Niveau 2->klev |
|---|
| 340 | DO k = 2, klev |
|---|
| 341 | DO i = 1, knon |
|---|
| 342 | flx_u_new(i,k) = Kcoefm(i,k)/RG/dtime * & |
|---|
| 343 | (u_new(i,k)-u_new(i,k-1)) |
|---|
| 344 | |
|---|
| 345 | flx_v_new(i,k) = Kcoefm(i,k)/RG/dtime * & |
|---|
| 346 | (v_new(i,k)-v_new(i,k-1)) |
|---|
| 347 | END DO |
|---|
| 348 | END DO |
|---|
| 349 | |
|---|
| 350 | !**************************************************************************************** |
|---|
| 351 | ! Calcul tendances |
|---|
| 352 | ! |
|---|
| 353 | !**************************************************************************************** |
|---|
| 354 | d_u_new(:,:) = 0.0 |
|---|
| 355 | d_v_new(:,:) = 0.0 |
|---|
| 356 | DO k = 1, klev |
|---|
| 357 | DO i = 1, knon |
|---|
| 358 | d_u_new(i,k) = u_new(i,k) - u_old(i,k) |
|---|
| 359 | d_v_new(i,k) = v_new(i,k) - v_old(i,k) |
|---|
| 360 | END DO |
|---|
| 361 | END DO |
|---|
| 362 | |
|---|
| 363 | END SUBROUTINE climb_wind_up |
|---|
| 364 | ! |
|---|
| 365 | !**************************************************************************************** |
|---|
| 366 | ! |
|---|
| 367 | END MODULE climb_wind_mod |
|---|
