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
Line 
1!> @file ls_forcing_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
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.
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!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: ls_forcing_mod.f90 2105 2017-01-06 16:02:14Z suehring $
27!
28! 2104 2017-01-06 16:01:15Z knoop
29! Bugfix for approximation related flux input conversion
30!
31! 2071 2016-11-17 11:22:14Z maronga
32! Small bugfix (Resler)
33!
34! 2037 2016-10-26 11:15:40Z knoop
35! Anelastic approximation implemented
36!
37! 2000 2016-08-20 18:09:15Z knoop
38! Forced header and separation lines into 80 columns
39!
40! 1850 2016-04-08 13:29:27Z maronga
41! Module renamed
42!
43!
44! 1682 2015-10-07 23:56:08Z knoop
45! Code annotations made doxygen readable
46!
47! 1602 2015-06-22 07:50:56Z heinze
48! PA0370 changed to PA0363
49!
50! 1382 2014-04-30 12:15:41Z boeske
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, ...
55!
56! 1365 2014-04-22 15:03:56Z boeske
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
62!
63! 1353 2014-04-08 15:21:23Z heinze
64! REAL constants provided with KIND-attribute
65!
66! 1320 2014-03-20 08:40:49Z raasch
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
72!
73! 1318 2014-03-17 13:35:16Z raasch
74! module interfaces removed
75!
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!
80! 1276 2014-01-15 13:40:41Z heinze
81! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
82!
83! 1249 2013-11-06 10:45:47Z heinze
84! remove call of user module
85! reformatting
86!
87! 1241 2013-10-30 11:36:58Z heinze
88! Initial revision
89!
90! Description:
91! ------------
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.
95!--------------------------------------------------------------------------------!
96 MODULE ls_forcing_mod
97 
98
99    PRIVATE
100    PUBLIC init_ls_forcing, ls_forcing_surf, ls_forcing_vert, ls_advec
101    SAVE
102
103    INTERFACE ls_advec
104       MODULE PROCEDURE ls_advec
105       MODULE PROCEDURE ls_advec_ij
106    END INTERFACE ls_advec
107
108 CONTAINS
109
110!------------------------------------------------------------------------------!
111! Description:
112! ------------
113!> @todo Missing subroutine description.
114!------------------------------------------------------------------------------!
115    SUBROUTINE init_ls_forcing
116
117       USE arrays_3d,                                                          &
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
121
122       USE control_parameters,                                                 &
123           ONLY:  end_time, lsf_surf, lsf_vert, message_string, nlsf
124
125       USE indices,                                                            &
126           ONLY:  ngp_sums_ls, nzb, nz, nzt
127
128       USE kinds
129
130       USE statistics,                                                         &
131           ONLY:  sums_ls_l
132
133
134       IMPLICIT NONE
135
136       CHARACTER(100) ::  chmess      !<
137       CHARACTER(1)   ::  hash        !<
138
139       INTEGER(iwp) ::  ierrn         !<
140       INTEGER(iwp) ::  finput = 90   !<
141       INTEGER(iwp) ::  k             !<
142       INTEGER(iwp) ::  nt             !<
143
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           !<
162
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),     &
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) )
170
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
176
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
182
183       ngp_sums_ls = (nz+2)*6
184
185       OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
186              FORM='FORMATTED', IOSTAT=ierrn )
187
188       IF ( ierrn /= 0 )  THEN
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
196       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
197       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
198       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
199
200       IF ( ierrn /= 0 )  THEN
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
207       nt     = 0
208       ierrn = 0
209
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)
215
216          IF ( ierrn < 0 )  THEN
217            WRITE ( message_string, * ) 'No time dependent surface variables ',&
218                              'in&LSF_DATA for end of run found'
219
220             CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
221          ENDIF
222       ENDDO
223
224       IF ( time_surf(1) > end_time )  THEN
225          WRITE ( message_string, * ) 'No time dependent surface variables in ',&
226                                     '&LSF_DATA for end of run found - ',      &
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 )
235          READ ( finput, *, IOSTAT = ierrn ) r_dummy
236       ENDDO
237
238!
239!--    Profiles of ug, vg and w_subs are read in (large scale forcing)
240
241       nt = 0
242       DO WHILE ( time_vert(nt) < end_time )
243          nt = nt + 1
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 ) ) 
250             READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
251             IF ( ierrn < 0 )  THEN
252                WRITE( message_string, * ) 'No time dependent vertical profiles',&
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
258          IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
259
260          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,  &
261                                           lowwsubs_vert, low_td_lsa_lpt,      &
262                                           low_td_lsa_q, low_td_sub_lpt,       &
263                                           low_td_sub_q
264          IF ( ierrn /= 0 )  THEN
265             message_string = 'errors in file LSF_DATA'
266             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
267          ENDIF
268
269          READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,            &
270                                           highvg_vert, highwsubs_vert,        &
271                                           high_td_lsa_lpt, high_td_lsa_q,     &
272                                           high_td_sub_lpt, high_td_sub_q
273     
274          IF ( ierrn /= 0 )  THEN
275             message_string = 'errors in file LSF_DATA'
276             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
277          ENDIF
278
279
280          DO  k = nzb, nzt+1
281             IF ( highheight < zu(k) )  THEN
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
290
291                ierrn = 0
292                READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,      &
293                                                 highvg_vert, highwsubs_vert,  &
294                                                 high_td_lsa_lpt,              &
295                                                 high_td_lsa_q,                &
296                                                 high_td_sub_lpt, high_td_sub_q
297
298                IF ( ierrn /= 0 )  THEN
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 )
304                ENDIF
305
306             ENDIF
307
308!
309!--          Interpolation of prescribed profiles in space
310             fac = (highheight-zu(k))/(highheight - lowheight)
311
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
318
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
327
328          ENDDO
329
330       ENDDO 
331
332!
333!--    Large scale vertical velocity has to be zero at the surface
334       wsubs_vert(nzb,:) = 0.0_wp
335 
336       IF ( time_vert(1) > end_time )  THEN
337          WRITE ( message_string, * ) 'Time dependent large scale profile ',   &
338                             'forcing from&LSF_DATA sets in after end of ' ,   &
339                             'simulation - lsf_vert is set to FALSE'
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
350!------------------------------------------------------------------------------!
351! Description:
352! ------------
353!> @todo Missing subroutine description.
354!------------------------------------------------------------------------------!
355    SUBROUTINE ls_forcing_surf ( time )
356
357       USE arrays_3d,                                                          &
358           ONLY:  p_surf, pt_surf, q_surf, qsws, qsws_surf, shf, shf_surf,     &
359                  heatflux_input_conversion, waterflux_input_conversion,       &
360                  time_surf, time_vert, ug, ug_vert, vg, vg_vert
361
362       USE control_parameters,                                                 &
363           ONLY:  bc_q_b, ibc_pt_b, ibc_q_b, pt_surface, q_surface,            &
364                  surface_pressure
365
366       USE indices,                                                            &
367           ONLY:  nzb
368
369       USE kinds
370
371       IMPLICIT NONE
372
373       INTEGER(iwp) ::  nt                     !<
374
375       REAL(wp)             :: fac            !<
376       REAL(wp), INTENT(in) :: time           !<
377
378!
379!--    Interpolation in time of LSF_DATA at the surface
380       nt = 1
381       DO WHILE ( time > time_surf(nt) )
382          nt = nt + 1
383       ENDDO
384       IF ( time /= time_surf(nt) )  THEN
385         nt = nt - 1
386       ENDIF
387
388       fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
389
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
394          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
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
400          shf        = ( shf_surf(nt) + fac * ( shf_surf(nt+1) - shf_surf(nt) )&
401                       ) * heatflux_input_conversion(nzb)
402          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
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
410          q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
411
412       ELSEIF ( ibc_q_b == 1 )  THEN
413
414          qsws = ( qsws_surf(nt) + fac * ( qsws_surf(nt+1) - qsws_surf(nt) )   &
415                 ) * waterflux_input_conversion(nzb)
416
417       ENDIF
418
419       surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
420
421    END SUBROUTINE ls_forcing_surf 
422
423
424!------------------------------------------------------------------------------!
425! Description:
426! ------------
427!> @todo Missing subroutine description.
428!------------------------------------------------------------------------------!
429    SUBROUTINE ls_forcing_vert ( time )
430
431       USE arrays_3d,                                                          &
432           ONLY:  time_vert, ug, ug_vert, vg, vg_vert, w_subs, wsubs_vert
433
434       USE control_parameters,                                                 &
435           ONLY:  large_scale_subsidence
436
437       USE kinds
438
439       IMPLICIT NONE
440
441       INTEGER(iwp) ::  nt                     !<
442
443       REAL(wp)             ::  fac           !<
444       REAL(wp), INTENT(in) ::  time          !<
445
446!
447!--    Interpolation in time of LSF_DATA for ug, vg and w_subs
448       nt = 1
449       DO WHILE ( time > time_vert(nt) )
450          nt = nt + 1
451       ENDDO
452       IF ( time /= time_vert(nt) )  THEN
453         nt = nt - 1
454       ENDIF
455
456       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
457
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) )
460
461       IF ( large_scale_subsidence )  THEN
462          w_subs = wsubs_vert(:,nt)                                            &
463                   + fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
464       ENDIF
465
466    END SUBROUTINE ls_forcing_vert
467
468
469!------------------------------------------------------------------------------!
470! Description:
471! ------------
472!> Call for all grid points
473!------------------------------------------------------------------------------!
474    SUBROUTINE ls_advec ( time, prog_var )
475
476       USE arrays_3d,                                                          &
477           ONLY:  td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, tend, time_vert       
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
488       CHARACTER (LEN=*) ::  prog_var   !<
489
490       REAL(wp), INTENT(in)  :: time    !<
491       REAL(wp) :: fac                  !< 
492
493       INTEGER(iwp) ::  i               !<
494       INTEGER(iwp) ::  j               !<
495       INTEGER(iwp) ::  k               !<
496       INTEGER(iwp) ::  nt               !<
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
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) )
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
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) )
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
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) )
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
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) )
562                      ENDDO
563                   ENDDO
564                ENDDO
565
566          END SELECT
567
568       ENDIF
569
570    END SUBROUTINE ls_advec
571
572
573!------------------------------------------------------------------------------!
574! Description:
575! ------------
576!> Call for grid point i,j
577!------------------------------------------------------------------------------!
578    SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
579
580       USE arrays_3d,                                                          &
581           ONLY:  td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, tend, time_vert       
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
592       CHARACTER (LEN=*) ::  prog_var   !<
593
594       REAL(wp), INTENT(in)  :: time    !<
595       REAL(wp) :: fac                  !<
596
597       INTEGER(iwp) ::  i               !<
598       INTEGER(iwp) ::  j               !<
599       INTEGER(iwp) ::  k               !<
600       INTEGER(iwp) ::  nt               !<
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
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) )
623             ENDDO
624
625          CASE ( 'q' )
626
627             DO  k = nzb_u_inner(j,i)+1, nzt
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) )
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
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) )
645                ENDDO
646 
647             CASE ( 'q' )
648
649                DO  k = nzb_u_inner(j,i)+1, nzt
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) )
652                ENDDO
653
654          END SELECT
655
656       ENDIF
657
658    END SUBROUTINE ls_advec_ij
659
660
661 END MODULE ls_forcing_mod
Note: See TracBrowser for help on using the repository browser.