source: palm/trunk/SOURCE/diffusion_u.f90 @ 4313

Last change on this file since 4313 was 4182, checked in by scharf, 22 months ago
  • corrected "Former revisions" section
  • minor formatting in "Former revisions" section
  • added "Author" section
  • Property svn:keywords set to Id
File size: 24.1 KB
RevLine 
[1873]1!> @file diffusion_u.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
[1036]9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
[3655]17! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1]21! -----------------
[1341]22!
[2233]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: diffusion_u.f90 4182 2019-08-22 15:20:23Z suehring $
[4182]27! Corrected "Former revisions" section
28!
29! 3655 2019-01-07 16:51:22Z knoop
[3634]30! OpenACC port for SPEC
[2716]31!
[4182]32! Revision 1.1  1997/09/12 06:23:51  raasch
33! Initial revision
34!
35!
[1]36! Description:
37! ------------
[1682]38!> Diffusion term of the u-component
39!> @todo additional damping (needed for non-cyclic bc) causes bad vectorization
40!>       and slows down the speed on NEC about 5-10%
[1]41!------------------------------------------------------------------------------!
[1682]42 MODULE diffusion_u_mod
43 
[1]44
45    PRIVATE
[2118]46    PUBLIC diffusion_u
[1]47
48    INTERFACE diffusion_u
49       MODULE PROCEDURE diffusion_u
50       MODULE PROCEDURE diffusion_u_ij
51    END INTERFACE diffusion_u
52
53 CONTAINS
54
55
56!------------------------------------------------------------------------------!
[1682]57! Description:
58! ------------
59!> Call for all grid points
[1]60!------------------------------------------------------------------------------!
[1001]61    SUBROUTINE diffusion_u
[1]62
[1320]63       USE arrays_3d,                                                          &
[2232]64           ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
[1320]65       
66       USE control_parameters,                                                 &
[2232]67           ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
[1320]68                  use_top_fluxes
69       
70       USE grid_variables,                                                     &
[2232]71           ONLY:  ddx, ddx2, ddy
[1320]72       
73       USE indices,                                                            &
[3241]74           ONLY:  nxlu, nxr, nyn, nys, nzb, nzt, wall_flags_0
[2232]75     
[1320]76       USE kinds
[1]77
[2232]78       USE surface_mod,                                                        &
79           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
80                   surf_usm_v
81
[1]82       IMPLICIT NONE
83
[2232]84       INTEGER(iwp) ::  i             !< running index x direction
85       INTEGER(iwp) ::  j             !< running index y direction
86       INTEGER(iwp) ::  k             !< running index z direction
87       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
88       INTEGER(iwp) ::  m             !< running index surface elements
[3547]89       INTEGER(iwp) ::  surf_e        !< end index of surface elements at (j,i)-gridpoint
90       INTEGER(iwp) ::  surf_s        !< start index of surface elements at (j,i)-gridpoint
[1001]91
[2232]92       REAL(wp)     ::  flag          !< flag to mask topography grid points
[3547]93       REAL(wp)     ::  kmym          !< diffusion coefficient on southward side of the u-gridbox - interpolated onto xu-yv grid
94       REAL(wp)     ::  kmyp          !< diffusion coefficient on northward side of the u-gridbox - interpolated onto xu-yv grid
95       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
96       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
[2232]97       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface       
98       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point
99       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point
100       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
[1]101
[56]102
[2232]103
[3634]104       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
105       !$ACC PRIVATE(surf_e, surf_s, flag, kmym, kmyp, kmzm, kmzp) &
106       !$ACC PRIVATE(mask_bottom, mask_north, mask_south, mask_top) &
107       !$ACC PRESENT(wall_flags_0, km) &
108       !$ACC PRESENT(u, v, w) &
109       !$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
110       !$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:1)) &
111       !$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:1)) &
112       !$ACC PRESENT(surf_usm_h, surf_usm_v(0:1)) &
113       !$ACC PRESENT(tend)
[106]114       DO  i = nxlu, nxr
[1001]115          DO  j = nys, nyn
[1]116!
117!--          Compute horizontal diffusion
[2232]118             DO  k = nzb+1, nzt
[1]119!
[2232]120!--             Predetermine flag to mask topography and wall-bounded grid points.
121!--             It is sufficient to masked only north- and south-facing surfaces, which
122!--             need special treatment for the u-component.
123                flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i),   1 ) ) 
124                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 1 ) )
125                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 1 ) )
126!
[1]127!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]128                kmyp = 0.25_wp *                                               &
[978]129                       ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
[1340]130                kmym = 0.25_wp *                                               &
[978]131                       ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
[1]132
[1320]133                tend(k,j,i) = tend(k,j,i)                                      &
[2232]134                        + 2.0_wp * (                                           &
135                                  km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )    &
136                                - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )    &
137                                   ) * ddx2 * flag                             &
138                        +          ( mask_north * (                            &
139                            kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy         &
140                          + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx         &
141                                                  )                            &
142                                   - mask_south * (                            &
143                            kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
144                          + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
145                                                  )                            &
146                                   ) * ddy  * flag                             
[1]147             ENDDO
148!
[2232]149!--          Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1)
150!--          surfaces. Note, in the the flat case, loops won't be entered as
151!--          start_index > end_index. Furtermore, note, no vertical natural surfaces
152!--          so far.           
153!--          Default-type surfaces
154             DO  l = 0, 1
155                surf_s = surf_def_v(l)%start_index(j,i)
156                surf_e = surf_def_v(l)%end_index(j,i)
157                DO  m = surf_s, surf_e
158                   k           = surf_def_v(l)%k(m)
159                   tend(k,j,i) = tend(k,j,i) +                                 &                   
160                                    surf_def_v(l)%mom_flux_uv(m) * ddy
161                ENDDO   
162             ENDDO
163!
164!--          Natural-type surfaces
165             DO  l = 0, 1
166                surf_s = surf_lsm_v(l)%start_index(j,i)
167                surf_e = surf_lsm_v(l)%end_index(j,i)
168                DO  m = surf_s, surf_e
169                   k           = surf_lsm_v(l)%k(m)
170                   tend(k,j,i) = tend(k,j,i) +                                 &                   
171                                    surf_lsm_v(l)%mom_flux_uv(m) * ddy
172                ENDDO   
173             ENDDO
174!
175!--          Urban-type surfaces
176             DO  l = 0, 1
177                surf_s = surf_usm_v(l)%start_index(j,i)
178                surf_e = surf_usm_v(l)%end_index(j,i)
179                DO  m = surf_s, surf_e
180                   k           = surf_usm_v(l)%k(m)
181                   tend(k,j,i) = tend(k,j,i) +                                 &                   
182                                    surf_usm_v(l)%mom_flux_uv(m) * ddy
183                ENDDO   
184             ENDDO
[51]185
[1]186!
[2232]187!--          Compute vertical diffusion. In case of simulating a surface layer,
188!--          respective grid diffusive fluxes are masked (flag 8) within this
189!--          loop, and added further below, else, simple gradient approach is
190!--          applied. Model top is also mask if top-momentum flux is given.
191             DO  k = nzb+1, nzt
[1]192!
[2232]193!--             Determine flags to mask topography below and above. Flag 1 is
194!--             used to mask topography in general, and flag 8 implies
195!--             information about use_surface_fluxes. Flag 9 is used to control
196!--             momentum flux at model top. 
197                mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
198                                     BTEST( wall_flags_0(k-1,j,i), 8 ) ) 
199                mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
200                                     BTEST( wall_flags_0(k+1,j,i), 8 ) ) *     &
201                              MERGE( 1.0_wp, 0.0_wp,                           &
202                                     BTEST( wall_flags_0(k+1,j,i), 9 ) ) 
203                flag        = MERGE( 1.0_wp, 0.0_wp,                           &
204                                     BTEST( wall_flags_0(k,j,i), 1 ) ) 
205!
[1]206!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]207                kmzp = 0.25_wp *                                               &
[1]208                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
[1340]209                kmzm = 0.25_wp *                                               &
[1]210                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
211
[1320]212                tend(k,j,i) = tend(k,j,i)                                      &
[2232]213                        + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
214                                   + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
215                                   ) * rho_air_zw(k)   * mask_top              &
216                          - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
217                                   + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
218                                   ) * rho_air_zw(k-1) * mask_bottom           &
219                          ) * ddzw(k) * drho_air(k) * flag
[1]220             ENDDO
221
222!
223!--          Vertical diffusion at the first grid point above the surface,
224!--          if the momentum flux at the bottom is given by the Prandtl law or
225!--          if it is prescribed by the user.
226!--          Difference quotient of the momentum flux is not formed over half
227!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
[1320]228!--          with other (LES) models showed that the values of the momentum
[1]229!--          flux becomes too large in this case.
230!--          The term containing w(k-1,..) (see above equation) is removed here
231!--          because the vertical velocity is assumed to be zero at the surface.
232             IF ( use_surface_fluxes )  THEN
233!
[2232]234!--             Default-type surfaces, upward-facing
235                surf_s = surf_def_h(0)%start_index(j,i)
236                surf_e = surf_def_h(0)%end_index(j,i)
237                DO  m = surf_s, surf_e
[1]238
[2232]239                   k   = surf_def_h(0)%k(m)
[1]240
[2232]241                   tend(k,j,i) = tend(k,j,i)                                   &
242                        + ( - ( - surf_def_h(0)%usws(m) )                      &
243                          ) * ddzw(k) * drho_air(k)
244                ENDDO
[102]245!
[2232]246!--             Default-type surfaces, dowward-facing
247                surf_s = surf_def_h(1)%start_index(j,i)
248                surf_e = surf_def_h(1)%end_index(j,i)
249                DO  m = surf_s, surf_e
250
251                   k   = surf_def_h(1)%k(m)
252
253                   tend(k,j,i) = tend(k,j,i)                                   &
254                        + ( - surf_def_h(1)%usws(m)                            &
255                          ) * ddzw(k) * drho_air(k)
256                ENDDO
[102]257!
[2232]258!--             Natural-type surfaces, upward-facing
259                surf_s = surf_lsm_h%start_index(j,i)
260                surf_e = surf_lsm_h%end_index(j,i)
261                DO  m = surf_s, surf_e
[102]262
[2232]263                   k   = surf_lsm_h%k(m)
264
265                   tend(k,j,i) = tend(k,j,i)                                   &
266                        + ( - ( - surf_lsm_h%usws(m) )                         &
267                          ) * ddzw(k) * drho_air(k)
268                ENDDO
269!
270!--             Urban-type surfaces, upward-facing
271                surf_s = surf_usm_h%start_index(j,i)
272                surf_e = surf_usm_h%end_index(j,i)
273                DO  m = surf_s, surf_e
274
275                   k   = surf_usm_h%k(m)
276
277                   tend(k,j,i) = tend(k,j,i)                                   &
278                       + ( - ( - surf_usm_h%usws(m) )                          &
279                         ) * ddzw(k) * drho_air(k)
280                ENDDO
281
[102]282             ENDIF
[2232]283!
284!--          Add momentum flux at model top
[2638]285             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[2232]286                surf_s = surf_def_h(2)%start_index(j,i)
287                surf_e = surf_def_h(2)%end_index(j,i)
288                DO  m = surf_s, surf_e
[102]289
[2232]290                   k   = surf_def_h(2)%k(m)
291
292                   tend(k,j,i) = tend(k,j,i)                                   &
293                        + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
294                ENDDO
295             ENDIF
296
[1]297          ENDDO
298       ENDDO
299
300    END SUBROUTINE diffusion_u
301
302
303!------------------------------------------------------------------------------!
[1682]304! Description:
305! ------------
306!> Call for grid point i,j
[1]307!------------------------------------------------------------------------------!
[1001]308    SUBROUTINE diffusion_u_ij( i, j )
[1]309
[1320]310       USE arrays_3d,                                                          &
[2232]311           ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
[1320]312       
313       USE control_parameters,                                                 &
[2232]314           ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
315                  use_top_fluxes
[1320]316       
317       USE grid_variables,                                                     &
[2232]318           ONLY:  ddx, ddx2, ddy
[1320]319       
320       USE indices,                                                            &
[2232]321           ONLY:  nzb, nzt, wall_flags_0
322     
[1320]323       USE kinds
[1]324
[2232]325       USE surface_mod,                                                        &
326           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
327                   surf_usm_v
328
[1]329       IMPLICIT NONE
330
[2232]331       INTEGER(iwp) ::  i             !< running index x direction
332       INTEGER(iwp) ::  j             !< running index y direction
333       INTEGER(iwp) ::  k             !< running index z direction
334       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
335       INTEGER(iwp) ::  m             !< running index surface elements
336       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
337       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
[1]338
[2232]339       REAL(wp)     ::  flag          !< flag to mask topography grid points
[3547]340       REAL(wp)     ::  kmym          !< diffusion coefficient on southward side of the u-gridbox - interpolated onto xu-yv grid
341       REAL(wp)     ::  kmyp          !<diffusion coefficient on northward side of the u-gridbox - interpolated onto xu-yv grid
342       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
343       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
[2232]344       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface       
345       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point
346       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point
347       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
348!
[1]349!--    Compute horizontal diffusion
[2232]350       DO  k = nzb+1, nzt
[1]351!
[2232]352!--       Predetermine flag to mask topography and wall-bounded grid points.
353!--       It is sufficient to masked only north- and south-facing surfaces, which
354!--       need special treatment for the u-component.
355          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i),   1 ) ) 
356          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 1 ) )
357          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 1 ) )
358!
[1]359!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]360          kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
361          kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
[1]362
[1320]363          tend(k,j,i) = tend(k,j,i)                                            &
[2232]364                       + 2.0_wp * (                                            &
365                                 km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )     &
366                               - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )     &
367                                   ) * ddx2 * flag                             &
368                                 + (                                           &
369                  mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i)     ) * ddy    &
370                                      + ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx    &
371                                      )                                        &
372                - mask_south * kmym * ( ( u(k,j,i)   - u(k,j-1,i)   ) * ddy    &
373                                      + ( v(k,j,i)   - v(k,j,i-1)   ) * ddx    &
374                                      )                                        &
375                                   ) * ddy  * flag
[1]376       ENDDO
377
378!
[2232]379!--    Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1)
380!--    surfaces. Note, in the the flat case, loops won't be entered as
381!--    start_index > end_index. Furtermore, note, no vertical natural surfaces
382!--    so far.           
383!--    Default-type surfaces
384       DO  l = 0, 1
385          surf_s = surf_def_v(l)%start_index(j,i)
386          surf_e = surf_def_v(l)%end_index(j,i)
387          DO  m = surf_s, surf_e
388             k           = surf_def_v(l)%k(m)
389             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
390          ENDDO   
391       ENDDO
[51]392!
[2232]393!--    Natural-type surfaces
394       DO  l = 0, 1
395          surf_s = surf_lsm_v(l)%start_index(j,i)
396          surf_e = surf_lsm_v(l)%end_index(j,i)
397          DO  m = surf_s, surf_e
398             k           = surf_lsm_v(l)%k(m)
399             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
400          ENDDO   
401       ENDDO
[1]402!
[2232]403!--    Urban-type surfaces
404       DO  l = 0, 1
405          surf_s = surf_usm_v(l)%start_index(j,i)
406          surf_e = surf_usm_v(l)%end_index(j,i)
407          DO  m = surf_s, surf_e
408             k           = surf_usm_v(l)%k(m)
409             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
410          ENDDO   
411       ENDDO
[1]412!
[2232]413!--    Compute vertical diffusion. In case of simulating a surface layer,
414!--    respective grid diffusive fluxes are masked (flag 8) within this
415!--    loop, and added further below, else, simple gradient approach is
416!--    applied. Model top is also mask if top-momentum flux is given.
417       DO  k = nzb+1, nzt
418!
419!--       Determine flags to mask topography below and above. Flag 1 is
420!--       used to mask topography in general, and flag 8 implies
421!--       information about use_surface_fluxes. Flag 9 is used to control
422!--       momentum flux at model top.
423          mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
424                               BTEST( wall_flags_0(k-1,j,i), 8 ) ) 
425          mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
426                               BTEST( wall_flags_0(k+1,j,i), 8 ) ) *           &
427                        MERGE( 1.0_wp, 0.0_wp,                                 &
428                               BTEST( wall_flags_0(k+1,j,i), 9 ) )
429          flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
430                               BTEST( wall_flags_0(k,j,i), 1 ) ) 
431!
[1]432!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]433          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
434          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
[1]435
[1320]436          tend(k,j,i) = tend(k,j,i)                                            &
[2232]437                        + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
438                                   + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
439                                   ) * rho_air_zw(k)   * mask_top              &
440                          - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
441                                   + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
442                                   ) * rho_air_zw(k-1) * mask_bottom           &
443                          ) * ddzw(k) * drho_air(k) * flag
[1]444       ENDDO
445
446!
[2232]447!--    Vertical diffusion at the first surface grid points, if the
[1]448!--    momentum flux at the bottom is given by the Prandtl law or if it is
449!--    prescribed by the user.
450!--    Difference quotient of the momentum flux is not formed over half of
451!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
[1320]452!--    other (LES) models showed that the values of the momentum flux becomes
[1]453!--    too large in this case.
454       IF ( use_surface_fluxes )  THEN
455!
[2232]456!--       Default-type surfaces, upward-facing
457          surf_s = surf_def_h(0)%start_index(j,i)
458          surf_e = surf_def_h(0)%end_index(j,i)
459          DO  m = surf_s, surf_e
[1]460
[2232]461             k   = surf_def_h(0)%k(m)
[1]462
[2232]463             tend(k,j,i) = tend(k,j,i)                                         &
464                        + ( - ( - surf_def_h(0)%usws(m) )                      &
465                          ) * ddzw(k) * drho_air(k)
466          ENDDO
[102]467!
[2232]468!--       Default-type surfaces, dowward-facing (except for model-top fluxes)
469          surf_s = surf_def_h(1)%start_index(j,i)
470          surf_e = surf_def_h(1)%end_index(j,i)
471          DO  m = surf_s, surf_e
472
473             k   = surf_def_h(1)%k(m)
474
475             tend(k,j,i) = tend(k,j,i)                                         &
476                        + ( - surf_def_h(1)%usws(m)                            &
477                          ) * ddzw(k) * drho_air(k)
478          ENDDO
[102]479!
[2232]480!--       Natural-type surfaces, upward-facing
481          surf_s = surf_lsm_h%start_index(j,i)
482          surf_e = surf_lsm_h%end_index(j,i)
483          DO  m = surf_s, surf_e
[102]484
[2232]485             k   = surf_lsm_h%k(m)
486
487             tend(k,j,i) = tend(k,j,i)                                         &
488                        + ( - ( - surf_lsm_h%usws(m) )                         &
489                          ) * ddzw(k) * drho_air(k)
490          ENDDO
491!
492!--       Urban-type surfaces, upward-facing
493          surf_s = surf_usm_h%start_index(j,i)
494          surf_e = surf_usm_h%end_index(j,i)
495          DO  m = surf_s, surf_e
496
497             k   = surf_usm_h%k(m)
498
499             tend(k,j,i) = tend(k,j,i)                                         &
500                       + ( - ( - surf_usm_h%usws(m) )                          &
501                         ) * ddzw(k) * drho_air(k)
502          ENDDO
503
[102]504       ENDIF
[2232]505!
506!--    Add momentum flux at model top
[2638]507       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[2232]508          surf_s = surf_def_h(2)%start_index(j,i)
509          surf_e = surf_def_h(2)%end_index(j,i)
510          DO  m = surf_s, surf_e
[102]511
[2232]512             k   = surf_def_h(2)%k(m)
513
514             tend(k,j,i) = tend(k,j,i)                                         &
515                           + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
516          ENDDO
517       ENDIF
518
519
[1]520    END SUBROUTINE diffusion_u_ij
521
522 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.