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

Last change on this file since 1365 was 1365, checked in by boeske, 8 years ago

large scale forcing enabled

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