source: palm/trunk/SOURCE/ls_forcing.f90 @ 1777

Last change on this file since 1777 was 1683, checked in by knoop, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 21.7 KB
RevLine 
[1682]1!> @file ls_forcing.f90
[1239]2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
[1310]16! Copyright 1997-2014 Leibniz Universitaet Hannover
[1239]17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
[1383]21!
[1683]22!
[1383]23! Former revisions:
24! -----------------
25! $Id: ls_forcing.f90 1683 2015-10-07 23:57:51Z hoffmann $
26!
[1683]27! 1682 2015-10-07 23:56:08Z knoop
28! Code annotations made doxygen readable
29!
[1603]30! 1602 2015-06-22 07:50:56Z heinze
31! PA0370 changed to PA0363
32!
[1383]33! 1382 2014-04-30 12:15:41Z boeske
[1382]34! Renamed variables which store large scale forcing tendencies
35! pt_lsa -> td_lsa_lpt, pt_subs -> td_sub_lpt,
36! q_lsa  -> td_lsa_q,   q_subs  -> td_sub_q,
37! high|lowpt_lsa -> high|low_td_lsa_lpt, ...
[1366]38!
39! 1365 2014-04-22 15:03:56Z boeske
[1365]40! Usage of large scale forcing for pt and q enabled:
41! Added new subroutine ls_advec for horizontal large scale advection and large
42! scale subsidence,
43! error message in init_ls_forcing specified,
44! variable t renamed nt
[1354]45!
46! 1353 2014-04-08 15:21:23Z heinze
47! REAL constants provided with KIND-attribute
48!
[1321]49! 1320 2014-03-20 08:40:49Z raasch
[1320]50! ONLY-attribute added to USE-statements,
51! kind-parameters added to all INTEGER and REAL declaration statements,
52! kinds are defined in new module kinds,
53! comment fields (!:) to be used for variable explanations added to
54! all variable declaration statements
[1239]55!
[1319]56! 1318 2014-03-17 13:35:16Z raasch
57! module interfaces removed
58!
[1300]59! 1299 2014-03-06 13:15:21Z heinze
60! Ensure a zero large scale vertical velocity at the surface
61! Bugfix: typo in case of boundary condition in if-clause
62!
[1277]63! 1276 2014-01-15 13:40:41Z heinze
64! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
65!
[1250]66! 1249 2013-11-06 10:45:47Z heinze
67! remove call of user module
68! reformatting
69!
[1242]70! 1241 2013-10-30 11:36:58Z heinze
71! Initial revision
[1239]72!
73! Description:
74! ------------
[1682]75!> Calculates large scale forcings (geostrophic wind and subsidence velocity) as
76!> well as surfaces fluxes dependent on time given in an external file (LSF_DATA).
77!> Code is based in parts on DALES and UCLA-LES.
[1239]78!--------------------------------------------------------------------------------!
[1682]79 MODULE ls_forcing_mod
80 
[1239]81
82    PRIVATE
[1365]83    PUBLIC init_ls_forcing, ls_forcing_surf, ls_forcing_vert, ls_advec
[1239]84    SAVE
85
[1365]86    INTERFACE ls_advec
87       MODULE PROCEDURE ls_advec
88       MODULE PROCEDURE ls_advec_ij
89    END INTERFACE ls_advec
[1239]90
91 CONTAINS
92
[1682]93!------------------------------------------------------------------------------!
94! Description:
95! ------------
96!> @todo Missing subroutine description.
97!------------------------------------------------------------------------------!
[1239]98    SUBROUTINE init_ls_forcing
99
[1320]100       USE arrays_3d,                                                          &
[1382]101           ONLY:  p_surf, pt_surf, q_surf, qsws_surf, shf_surf, td_lsa_lpt,    &
102                  td_lsa_q, td_sub_lpt, td_sub_q, time_surf, time_vert,        &
103                  ug_vert, vg_vert, wsubs_vert, zu
[1239]104
[1320]105       USE control_parameters,                                                 &
106           ONLY:  end_time, lsf_surf, lsf_vert, message_string, nlsf
107
108       USE indices,                                                            &
[1365]109           ONLY:  ngp_sums_ls, nzb, nz, nzt
[1320]110
111       USE kinds
112
[1365]113       USE statistics,                                                         &
114           ONLY:  sums_ls_l
[1320]115
[1365]116
[1239]117       IMPLICIT NONE
118
[1682]119       CHARACTER(100) ::  chmess      !<
120       CHARACTER(1)   ::  hash        !<
[1239]121
[1682]122       INTEGER(iwp) ::  ierrn         !<
123       INTEGER(iwp) ::  finput = 90   !<
124       INTEGER(iwp) ::  k             !<
125       INTEGER(iwp) ::  nt             !<
[1320]126
[1682]127       REAL(wp) ::  fac               !<
128       REAL(wp) ::  highheight        !<
129       REAL(wp) ::  highug_vert       !<
130       REAL(wp) ::  highvg_vert       !<
131       REAL(wp) ::  highwsubs_vert    !<
132       REAL(wp) ::  lowheight         !<
133       REAL(wp) ::  lowug_vert        !<
134       REAL(wp) ::  lowvg_vert        !<
135       REAL(wp) ::  lowwsubs_vert     !<
136       REAL(wp) ::  high_td_lsa_lpt   !<
137       REAL(wp) ::  low_td_lsa_lpt    !<
138       REAL(wp) ::  high_td_lsa_q     !<
139       REAL(wp) ::  low_td_lsa_q      !<
140       REAL(wp) ::  high_td_sub_lpt   !<
141       REAL(wp) ::  low_td_sub_lpt    !<
142       REAL(wp) ::  high_td_sub_q     !<
143       REAL(wp) ::  low_td_sub_q      !<
144       REAL(wp) ::  r_dummy           !<
[1320]145
[1382]146       ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),              &
147                 qsws_surf(0:nlsf), shf_surf(0:nlsf),                          &
148                 td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),     &
149                 td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),     &
[1365]150                 time_vert(0:nlsf), time_surf(0:nlsf),                         &
151                 ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),         &
152                 wsubs_vert(nzb:nzt+1,0:nlsf) )
[1239]153
[1382]154       p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
155       shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
156       td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
157       time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
158       wsubs_vert = 0.0_wp
[1239]159
[1365]160!
161!--    Array for storing large scale forcing and nudging tendencies at each
162!--    timestep for data output
163       ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
164       sums_ls_l = 0.0_wp
[1239]165
[1365]166       ngp_sums_ls = (nz+2)*6
167
[1249]168       OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
169              FORM='FORMATTED', IOSTAT=ierrn )
[1239]170
[1249]171       IF ( ierrn /= 0 )  THEN
[1239]172          message_string = 'file LSF_DATA does not exist'
173          CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
174       ENDIF
175
176       ierrn = 0
177!
178!--    First three lines of LSF_DATA contain header
[1249]179       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
180       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
181       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
[1239]182
[1249]183       IF ( ierrn /= 0 )  THEN
[1239]184          message_string = 'errors in file LSF_DATA'
185          CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
186       ENDIF
187
188!
189!--    Surface values are read in
[1365]190       nt     = 0
[1239]191       ierrn = 0
192
[1365]193       DO WHILE ( time_surf(nt) < end_time )
194          nt = nt + 1
195          READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),      &
196                                             qsws_surf(nt), pt_surf(nt),       &
197                                             q_surf(nt), p_surf(nt)
[1239]198
[1249]199          IF ( ierrn < 0 )  THEN
[1299]200            WRITE ( message_string, * ) 'No time dependent surface variables ',&
[1239]201                              'in&LSF_DATA for end of run found'
202
[1602]203             CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
[1239]204          ENDIF
205       ENDDO
206
207
[1249]208       IF ( time_surf(1) > end_time )  THEN
[1299]209          WRITE ( message_string, * ) 'No time dependent surface variables in ',&
[1365]210                                     '&LSF_DATA for end of run found - ',      &
[1239]211                                     'lsf_surf is set to FALSE'
212          CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
213          lsf_surf = .FALSE.
214       ENDIF
215
216!
217!--    Go to the end of the list with surface variables
218       DO WHILE ( ierrn == 0 )
[1249]219          READ ( finput, *, IOSTAT = ierrn ) r_dummy
[1239]220       ENDDO
221
222!
223!--    Profiles of ug, vg and w_subs are read in (large scale forcing)
224
[1365]225       nt = 0
226       DO WHILE ( time_vert(nt) < end_time )
227          nt = nt + 1
[1239]228          hash = "#"
229          ierrn = 1 ! not zero
230!
231!--       Search for the next line consisting of "# time",
232!--       from there onwards the profiles will be read
233          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
[1365]234             READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
[1249]235             IF ( ierrn < 0 )  THEN
[1299]236                WRITE( message_string, * ) 'No time dependent vertical profiles',&
[1239]237                                 ' in&LSF_DATA for end of run found'
238                CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
239             ENDIF
240          ENDDO
241
[1365]242          IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
[1239]243
[1365]244          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,  &
[1382]245                                           lowwsubs_vert, low_td_lsa_lpt,      &
246                                           low_td_lsa_q, low_td_sub_lpt,       &
247                                           low_td_sub_q
[1249]248          IF ( ierrn /= 0 )  THEN
[1239]249             message_string = 'errors in file LSF_DATA'
[1299]250             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
[1239]251          ENDIF
252
[1365]253          READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,            &
254                                           highvg_vert, highwsubs_vert,        &
[1382]255                                           high_td_lsa_lpt, high_td_lsa_q,     &
256                                           high_td_sub_lpt, high_td_sub_q
[1239]257     
[1249]258          IF ( ierrn /= 0 )  THEN
[1239]259             message_string = 'errors in file LSF_DATA'
[1299]260             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
[1239]261          ENDIF
262
[1241]263
[1239]264          DO  k = nzb, nzt+1
[1249]265             IF ( highheight < zu(k) )  THEN
[1382]266                lowheight      = highheight
267                lowug_vert     = highug_vert
268                lowvg_vert     = highvg_vert
269                lowwsubs_vert  = highwsubs_vert
270                low_td_lsa_lpt = high_td_lsa_lpt
271                low_td_lsa_q   = high_td_lsa_q
272                low_td_sub_lpt = high_td_sub_lpt
273                low_td_sub_q   = high_td_sub_q
[1239]274
275                ierrn = 0
[1365]276                READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,      &
277                                                 highvg_vert, highwsubs_vert,  &
[1382]278                                                 high_td_lsa_lpt,              &
279                                                 high_td_lsa_q,                &
280                                                 high_td_sub_lpt, high_td_sub_q
[1239]281
[1249]282                IF ( ierrn /= 0 )  THEN
[1365]283                   WRITE( message_string, * ) 'zu(nzt+1) = ', zu(nzt+1), 'm ', &
284                        'is higher than the maximum height in LSF_DATA which ',&
285                        'is ', lowheight, 'm. Interpolation on PALM ',         &
286                        'grid is not possible.'
287                   CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
[1239]288                ENDIF
289
290             ENDIF
291
292!
293!--          Interpolation of prescribed profiles in space
294             fac = (highheight-zu(k))/(highheight - lowheight)
295
[1382]296             ug_vert(k,nt)    = fac * lowug_vert                               &
297                                + ( 1.0_wp - fac ) * highug_vert
298             vg_vert(k,nt)    = fac * lowvg_vert                               &
299                                + ( 1.0_wp - fac ) * highvg_vert
300             wsubs_vert(k,nt) = fac * lowwsubs_vert                            &
301                                + ( 1.0_wp - fac ) * highwsubs_vert
[1239]302
[1382]303             td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                           &
304                                + ( 1.0_wp - fac ) * high_td_lsa_lpt
305             td_lsa_q(k,nt)   = fac * low_td_lsa_q                             &
306                                + ( 1.0_wp - fac ) * high_td_lsa_q
307             td_sub_lpt(k,nt) = fac * low_td_sub_lpt                           &
308                                + ( 1.0_wp - fac ) * high_td_sub_lpt
309             td_sub_q(k,nt)   = fac * low_td_sub_q                             &
310                                + ( 1.0_wp - fac ) * high_td_sub_q
[1365]311
[1239]312          ENDDO
313
314       ENDDO 
315
[1299]316!
317!--    Large scale vertical velocity has to be zero at the surface
[1353]318       wsubs_vert(nzb,:) = 0.0_wp
[1299]319 
[1249]320       IF ( time_vert(1) > end_time )  THEN
[1365]321          WRITE ( message_string, * ) 'Time dependent large scale profile ',   &
322                             'forcing from&LSF_DATA sets in after end of ' ,   &
[1249]323                             'simulation - lsf_vert is set to FALSE'
[1239]324          CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
325          lsf_vert = .FALSE.
326       ENDIF
327
328       CLOSE( finput )
329
330
331    END SUBROUTINE init_ls_forcing 
332
333
[1682]334!------------------------------------------------------------------------------!
335! Description:
336! ------------
337!> @todo Missing subroutine description.
338!------------------------------------------------------------------------------!
[1239]339    SUBROUTINE ls_forcing_surf ( time )
340
[1320]341       USE arrays_3d,                                                          &
342           ONLY:  p_surf, pt_surf, q_surf, qsws, qsws_surf, shf, shf_surf,     &
343                  time_surf, time_vert, ug, ug_vert, vg, vg_vert
[1239]344
[1320]345       USE control_parameters,                                                 &
346           ONLY:  bc_q_b, ibc_pt_b, ibc_q_b, pt_surface, q_surface,            &
347                  surface_pressure
348
349       USE kinds
350
[1239]351       IMPLICIT NONE
352
[1682]353       INTEGER(iwp) ::  nt                     !<
[1239]354
[1682]355       REAL(wp)             :: fac            !<
356       REAL(wp), INTENT(in) :: time           !<
[1320]357
[1239]358!
359!--    Interpolation in time of LSF_DATA at the surface
[1365]360       nt = 1
361       DO WHILE ( time > time_surf(nt) )
362          nt = nt + 1
[1239]363       ENDDO
[1365]364       IF ( time /= time_surf(nt) )  THEN
365         nt = nt - 1
[1239]366       ENDIF
367
[1365]368       fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
[1239]369
[1276]370       IF ( ibc_pt_b == 0 )  THEN
371!
372!--       In case of Dirichlet boundary condition shf must not
373!--       be set - it is calculated via MOST in prandtl_fluxes
[1365]374          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
[1276]375
376       ELSEIF ( ibc_pt_b == 1 )  THEN
377!
378!--       In case of Neumann boundary condition pt_surface is needed for
379!--       calculation of reference density
[1365]380          shf        = shf_surf(nt) + fac * ( shf_surf(nt+1) - shf_surf(nt) )
381          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
[1276]382
383       ENDIF
384
385       IF ( ibc_q_b == 0 )  THEN
386!
387!--       In case of Dirichlet boundary condition qsws must not
388!--       be set - it is calculated via MOST in prandtl_fluxes
[1365]389          q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
[1276]390
[1299]391       ELSEIF ( ibc_q_b == 1 )  THEN
[1276]392
[1365]393          qsws = qsws_surf(nt) + fac * ( qsws_surf(nt+1) - qsws_surf(nt) )
[1276]394
395       ENDIF
396
[1365]397       surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
[1239]398
399    END SUBROUTINE ls_forcing_surf 
400
401
[1682]402!------------------------------------------------------------------------------!
403! Description:
404! ------------
405!> @todo Missing subroutine description.
406!------------------------------------------------------------------------------!
[1239]407    SUBROUTINE ls_forcing_vert ( time )
408
[1320]409       USE arrays_3d,                                                          &
410           ONLY:  time_vert, ug, ug_vert, vg, vg_vert, w_subs, wsubs_vert
[1239]411
[1320]412       USE control_parameters,                                                 &
413           ONLY:  large_scale_subsidence
414
415       USE kinds
416
[1239]417       IMPLICIT NONE
418
[1682]419       INTEGER(iwp) ::  nt                     !<
[1239]420
[1682]421       REAL(wp)             ::  fac           !<
422       REAL(wp), INTENT(in) ::  time          !<
[1320]423
[1239]424!
425!--    Interpolation in time of LSF_DATA for ug, vg and w_subs
[1365]426       nt = 1
427       DO WHILE ( time > time_vert(nt) )
428          nt = nt + 1
[1239]429       ENDDO
[1365]430       IF ( time /= time_vert(nt) )  THEN
431         nt = nt - 1
[1239]432       ENDIF
433
[1365]434       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
[1239]435
[1365]436       ug     = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) )
437       vg     = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) )
[1239]438
[1249]439       IF ( large_scale_subsidence )  THEN
[1365]440          w_subs = wsubs_vert(:,nt)                                            &
441                   + fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
[1239]442       ENDIF
443
444    END SUBROUTINE ls_forcing_vert
445
446
[1365]447!------------------------------------------------------------------------------!
[1682]448! Description:
449! ------------
450!> Call for all grid points
[1365]451!------------------------------------------------------------------------------!
452    SUBROUTINE ls_advec ( time, prog_var )
453
454       USE arrays_3d,                                                          &
[1382]455           ONLY:  td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, tend, time_vert       
[1365]456       
457       USE control_parameters,                                                 &
458           ONLY:  large_scale_subsidence, use_subsidence_tendencies
459       
460       USE indices
461       
462       USE kinds
463
464       IMPLICIT NONE
465
[1682]466       CHARACTER (LEN=*) ::  prog_var   !<
[1365]467
[1682]468       REAL(wp), INTENT(in)  :: time    !<
469       REAL(wp) :: fac                  !< 
[1365]470
[1682]471       INTEGER(iwp) ::  i               !<
472       INTEGER(iwp) ::  j               !<
473       INTEGER(iwp) ::  k               !<
474       INTEGER(iwp) ::  nt               !<
[1365]475
476!
477!--    Interpolation in time of LSF_DATA
478       nt = 1
479       DO WHILE ( time > time_vert(nt) )
480          nt = nt + 1
481       ENDDO
482       IF ( time /= time_vert(nt) )  THEN
483         nt = nt - 1
484       ENDIF
485
486       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
487
488!
489!--    Add horizontal large scale advection tendencies of pt and q
490       SELECT CASE ( prog_var )
491
492          CASE ( 'pt' )
493
494             DO  i = nxl, nxr
495                DO  j = nys, nyn
496                   DO  k = nzb_u_inner(j,i)+1, nzt
[1382]497                      tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) + fac *     &
498                                    ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
[1365]499                   ENDDO
500                ENDDO
501             ENDDO
502
503          CASE ( 'q' )
504
505             DO  i = nxl, nxr
506                DO  j = nys, nyn
507                   DO  k = nzb_u_inner(j,i)+1, nzt
[1382]508                      tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) + fac *       &
509                                    ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
[1365]510                   ENDDO
511                ENDDO
512             ENDDO
513
514       END SELECT
515
516!
517!--    Subsidence of pt and q with prescribed subsidence tendencies
518       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
519
520          SELECT CASE ( prog_var )
521
522             CASE ( 'pt' )
523
524                DO  i = nxl, nxr
525                   DO  j = nys, nyn
526                      DO  k = nzb_u_inner(j,i)+1, nzt
[1382]527                         tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac *  &
528                                       ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
[1365]529                      ENDDO
530                   ENDDO
531                ENDDO
532 
533             CASE ( 'q' )
534
535                DO  i = nxl, nxr
536                   DO  j = nys, nyn
537                      DO  k = nzb_u_inner(j,i)+1, nzt
[1382]538                         tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac *    &
539                                       ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
[1365]540                      ENDDO
541                   ENDDO
542                ENDDO
543
544          END SELECT
545
546       ENDIF
547
548    END SUBROUTINE ls_advec
549
550
551!------------------------------------------------------------------------------!
[1682]552! Description:
553! ------------
554!> Call for grid point i,j
[1365]555!------------------------------------------------------------------------------!
556    SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
557
558       USE arrays_3d,                                                          &
[1382]559           ONLY:  td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, tend, time_vert       
[1365]560       
561       USE control_parameters,                                                 &
562           ONLY:  large_scale_subsidence, use_subsidence_tendencies
563       
564       USE indices
565       
566       USE kinds
567
568       IMPLICIT NONE
569
[1682]570       CHARACTER (LEN=*) ::  prog_var   !<
[1365]571
[1682]572       REAL(wp), INTENT(in)  :: time    !<
573       REAL(wp) :: fac                  !<
[1365]574
[1682]575       INTEGER(iwp) ::  i               !<
576       INTEGER(iwp) ::  j               !<
577       INTEGER(iwp) ::  k               !<
578       INTEGER(iwp) ::  nt               !<
[1365]579
580!
581!--    Interpolation in time of LSF_DATA
582       nt = 1
583       DO WHILE ( time > time_vert(nt) )
584          nt = nt + 1
585       ENDDO
586       IF ( time /= time_vert(nt) )  THEN
587         nt = nt - 1
588       ENDIF
589
590       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
591
592!
593!--    Add horizontal large scale advection tendencies of pt and q
594       SELECT CASE ( prog_var )
595
596          CASE ( 'pt' )
597
598             DO  k = nzb_u_inner(j,i)+1, nzt
[1382]599                tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt)                   &
600                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
[1365]601             ENDDO
602
603          CASE ( 'q' )
604
605             DO  k = nzb_u_inner(j,i)+1, nzt
[1382]606                tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt)                     &
607                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
[1365]608             ENDDO
609
610       END SELECT
611
612!
613!--    Subsidence of pt and q with prescribed profiles
614       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
615
616          SELECT CASE ( prog_var )
617
618             CASE ( 'pt' )
619
620                DO  k = nzb_u_inner(j,i)+1, nzt
[1382]621                   tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac *        &
622                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
[1365]623                ENDDO
624 
625             CASE ( 'q' )
626
627                DO  k = nzb_u_inner(j,i)+1, nzt
[1382]628                   tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac *          &
629                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
[1365]630                ENDDO
631
632          END SELECT
633
634       ENDIF
635
636    END SUBROUTINE ls_advec_ij
637
638
[1239]639 END MODULE ls_forcing_mod
Note: See TracBrowser for help on using the repository browser.