source: palm/trunk/SOURCE/ls_forcing_mod.f90 @ 2106

Last change on this file since 2106 was 2105, checked in by knoop, 8 years ago

last commit documented

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