source: palm/trunk/SOURCE/nudging_mod.f90 @ 1952

Last change on this file since 1952 was 1851, checked in by maronga, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 19.9 KB
Line 
1!> @file nudging_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 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!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: nudging_mod.f90 1851 2016-04-08 13:32:50Z suehring $
26!
27! 1850 2016-04-08 13:29:27Z maronga
28! Module renamed
29!
30!
31! 1757 2016-02-22 15:49:32Z maronga
32! Bugfix: allow for using higher vertical resolution in nudging file than grid
33! spacing in the LES model
34!
35! 1682 2015-10-07 23:56:08Z knoop
36! Code annotations made doxygen readable
37!
38! 1398 2014-05-07 11:15:00Z heinze
39! Subroutine nudge_ref is extended to set u_init and v_init to the current
40! nudging profiles
41!
42! 1382 2014-04-30 12:15:41Z boeske
43! Changed the weighting factor that is used in the summation of nudging
44! tendencies for profile data output from weight_pres to weight_substep,
45! added Neumann boundary conditions for profile data output of nudging terms at
46! nzt+1
47!
48! 1380 2014-04-28 12:40:45Z heinze
49! Subroutine nudge_ref added to account for proper upper scalar boundary
50! conditions in case of nudging
51!
52! 1365 2014-04-22 15:03:56Z boeske
53! Variable t renamed nt, variable currtnudge renamed tmp_tnudge,
54! summation of nudging tendencies for data output added
55! +sums_ls_l, tmp_tend
56! Added new subroutine calc_tnudge, which calculates the current nudging time
57! scale at each time step
58!
59! 1355 2014-04-10 10:21:29Z heinze
60! Error message specified.
61!
62! 1353 2014-04-08 15:21:23Z heinze
63! REAL constants provided with KIND-attribute
64!
65! 1320 2014-03-20 08:40:49Z raasch
66! ONLY-attribute added to USE-statements,
67! kind-parameters added to all INTEGER and REAL declaration statements,
68! kinds are defined in new module kinds,
69! old module precision_kind is removed,
70! revision history before 2012 removed,
71! comment fields (!:) to be used for variable explanations added to
72! all variable declaration statements
73!
74! 1318 2014-03-17 13:35:16Z raasch
75! module interfaces removed
76!
77! 1268 2013-12-12 09:47:53Z heinze
78! bugfix: argument of calc_mean_profile corrected
79!
80! 1251 2013-11-07 08:14:30Z heinze
81! bugfix: calculate dtm and dtp also in vector version
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!> Nudges u, v, pt and q to given profiles on a relaxation timescale tnudge.
93!> Profiles are read in from NUDGIN_DATA. Code is based on Neggers et al. (2012)
94!> and also part of DALES and UCLA-LES.
95!--------------------------------------------------------------------------------!
96 MODULE nudge_mod
97 
98
99    PRIVATE
100    PUBLIC init_nudge, calc_tnudge, nudge, nudge_ref
101    SAVE
102
103    INTERFACE nudge
104       MODULE PROCEDURE nudge
105       MODULE PROCEDURE nudge_ij
106    END INTERFACE nudge
107
108 CONTAINS
109
110!------------------------------------------------------------------------------!
111! Description:
112! ------------
113!> @todo Missing subroutine description.
114!------------------------------------------------------------------------------!
115    SUBROUTINE init_nudge
116
117       USE arrays_3d,                                                          &
118           ONLY:  ptnudge, qnudge, timenudge, tmp_tnudge, tnudge, unudge,      &
119                  vnudge, wnudge, zu
120
121       USE control_parameters,                                                 &
122           ONLY:  dt_3d, lptnudge, lqnudge, lunudge, lvnudge, lwnudge,         &
123                   message_string, ntnudge
124
125       USE indices,                                                            &
126           ONLY:  nzb, nzt
127
128       USE kinds
129
130       IMPLICIT NONE
131
132
133       INTEGER(iwp) ::  finput = 90  !<
134       INTEGER(iwp) ::  ierrn        !<
135       INTEGER(iwp) ::  k            !<
136       INTEGER(iwp) ::  nt            !<
137
138       CHARACTER(1) ::  hash     !<
139
140       REAL(wp) ::  highheight   !<
141       REAL(wp) ::  highqnudge   !<
142       REAL(wp) ::  highptnudge  !<
143       REAL(wp) ::  highunudge   !<
144       REAL(wp) ::  highvnudge   !<
145       REAL(wp) ::  highwnudge   !<
146       REAL(wp) ::  hightnudge   !<
147
148       REAL(wp) ::  lowheight    !<
149       REAL(wp) ::  lowqnudge    !<
150       REAL(wp) ::  lowptnudge   !<
151       REAL(wp) ::  lowunudge    !<
152       REAL(wp) ::  lowvnudge    !<
153       REAL(wp) ::  lowwnudge    !<
154       REAL(wp) ::  lowtnudge    !<
155
156       REAL(wp) ::  fac          !<
157
158       ALLOCATE( ptnudge(nzb:nzt+1,1:ntnudge), qnudge(nzb:nzt+1,1:ntnudge), &
159                 tnudge(nzb:nzt+1,1:ntnudge), unudge(nzb:nzt+1,1:ntnudge),  &
160                 vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge)  )
161
162       ALLOCATE( tmp_tnudge(nzb:nzt) )
163
164       ALLOCATE( timenudge(0:ntnudge) )
165
166       ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
167       vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
168!
169!--    Initialize array tmp_nudge with a current nudging time scale of 6 hours
170       tmp_tnudge = 21600.0_wp
171
172       nt = 0
173       OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
174              FORM='FORMATTED', IOSTAT=ierrn )
175
176       IF ( ierrn /= 0 )  THEN
177          message_string = 'file NUDGING_DATA does not exist'
178          CALL message( 'nudging', 'PA0365', 1, 2, 0, 6, 0 )
179       ENDIF
180
181       ierrn = 0
182
183 rloop:DO
184          nt = nt + 1
185          hash = "#"
186          ierrn = 1 ! not zero
187!
188!--       Search for the next line consisting of "# time",
189!--       from there onwards the profiles will be read
190          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
191         
192            READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
193            IF ( ierrn < 0 )  EXIT rloop
194
195          ENDDO
196
197          ierrn = 0
198          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowtnudge, lowunudge,   &
199                                           lowvnudge, lowwnudge , lowptnudge, &
200                                           lowqnudge
201
202          IF ( ierrn /= 0 )  THEN
203             message_string = 'errors in file NUDGING_DATA'
204             CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
205          ENDIF
206
207          ierrn = 0
208          READ ( finput, *, IOSTAT=ierrn ) highheight, hightnudge, highunudge,   &
209                                           highvnudge, highwnudge , highptnudge, &
210                                           highqnudge
211
212          IF ( ierrn /= 0 )  THEN
213             message_string = 'errors in file NUDGING_DATA'
214             CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
215          ENDIF
216
217          DO  k = nzb, nzt+1
218             DO WHILE ( highheight < zu(k) )
219                lowheight  = highheight
220                lowtnudge  = hightnudge
221                lowunudge  = highunudge
222                lowvnudge  = highvnudge
223                lowwnudge  = highwnudge
224                lowptnudge = highptnudge
225                lowqnudge  = highqnudge
226 
227                ierrn = 0
228                READ ( finput, *, IOSTAT=ierrn )  highheight , hightnudge , &
229                                                  highunudge , highvnudge , &
230                                                  highwnudge , highptnudge, &
231                                                  highqnudge
232                IF (ierrn /= 0 )  THEN
233                   WRITE( message_string, * ) 'zu(nzt+1) = ', zu(nzt+1), 'm is ',&
234                        'higher than the maximum height in NUDING_DATA which ',  &
235                        'is ', lowheight, 'm. Interpolation on PALM ',           &
236                        'grid is not possible.'
237                   CALL message( 'nudging', 'PA0364', 1, 2, 0, 6, 0 )
238                ENDIF
239             ENDDO
240
241!
242!--          Interpolation of prescribed profiles in space
243
244             fac = ( highheight - zu(k) ) / ( highheight - lowheight )
245
246             tnudge(k,nt)  = fac * lowtnudge  + ( 1.0_wp - fac ) * hightnudge
247             unudge(k,nt)  = fac * lowunudge  + ( 1.0_wp - fac ) * highunudge
248             vnudge(k,nt)  = fac * lowvnudge  + ( 1.0_wp - fac ) * highvnudge
249             wnudge(k,nt)  = fac * lowwnudge  + ( 1.0_wp - fac ) * highwnudge
250             ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
251             qnudge(k,nt)  = fac * lowqnudge  + ( 1.0_wp - fac ) * highqnudge
252          ENDDO
253
254       ENDDO rloop
255
256       CLOSE ( finput )
257
258!
259!--    Prevent nudging if nudging profiles exhibt too small values
260!--    not used so far
261       lptnudge  = ANY( ABS( ptnudge ) > 1.0e-8_wp )
262       lqnudge   = ANY( ABS( qnudge )  > 1.0e-8_wp )
263       lunudge   = ANY( ABS( unudge )  > 1.0e-8_wp )
264       lvnudge   = ANY( ABS( vnudge )  > 1.0e-8_wp )
265       lwnudge   = ANY( ABS( wnudge )  > 1.0e-8_wp )
266
267    END SUBROUTINE init_nudge
268
269
270!------------------------------------------------------------------------------!
271! Description:
272! ------------
273!> @todo Missing subroutine description.
274!------------------------------------------------------------------------------!
275    SUBROUTINE calc_tnudge ( time )
276
277       USE arrays_3d,                                                          &
278           ONLY:  timenudge, tmp_tnudge, tnudge
279
280       USE control_parameters,                                                 &
281           ONLY:  dt_3d 
282
283       USE indices,                                                            &
284           ONLY:  nzb, nzt
285
286       USE kinds
287
288       IMPLICIT NONE
289
290
291       REAL(wp) ::  dtm         !<
292       REAL(wp) ::  dtp         !<
293       REAL(wp) ::  time        !<
294
295       INTEGER(iwp) ::  k   !<
296       INTEGER(iwp) ::  nt  !<
297
298       nt = 1
299       DO WHILE ( time > timenudge(nt) )
300         nt = nt+1
301       ENDDO
302       IF ( time /= timenudge(1) ) THEN
303         nt = nt-1
304       ENDIF
305
306       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
307       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
308
309       DO  k = nzb, nzt
310          tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
311       ENDDO
312
313    END SUBROUTINE calc_tnudge
314
315!------------------------------------------------------------------------------!
316! Description:
317! ------------
318!> Call for all grid points
319!------------------------------------------------------------------------------!
320    SUBROUTINE nudge ( time, prog_var )
321
322       USE arrays_3d,                                                          &
323           ONLY:  pt, ptnudge, q, qnudge, tend, timenudge, tmp_tnudge, tnudge, &
324                  u, unudge, v, vnudge
325
326       USE control_parameters,                                                 &
327           ONLY:  dt_3d, intermediate_timestep_count, message_string
328
329       USE indices,                                                            &
330           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_u_inner, nzt
331
332       USE kinds
333
334       USE statistics,                                                         &
335           ONLY:  hom, sums_ls_l, weight_substep
336
337       IMPLICIT NONE
338
339       CHARACTER (LEN=*) ::  prog_var  !<
340
341       REAL(wp) ::  tmp_tend    !<
342       REAL(wp) ::  dtm         !<
343       REAL(wp) ::  dtp         !<
344       REAL(wp) ::  time        !<
345
346       INTEGER(iwp) ::  i  !<
347       INTEGER(iwp) ::  j  !<
348       INTEGER(iwp) ::  k  !<
349       INTEGER(iwp) ::  nt  !<
350
351
352       nt = 1
353       DO WHILE ( time > timenudge(nt) )
354         nt = nt+1
355       ENDDO
356       IF ( time /= timenudge(1) ) THEN
357         nt = nt-1
358       ENDIF
359
360       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
361       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
362
363       SELECT CASE ( prog_var )
364
365          CASE ( 'u' )
366
367             DO  i = nxl, nxr
368                DO  j = nys, nyn
369
370                   DO  k = nzb_u_inner(j,i)+1, nzt
371
372                      tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +     &
373                                     unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
374
375                      tend(k,j,i) = tend(k,j,i) + tmp_tend
376
377                      sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend *             &
378                                     weight_substep(intermediate_timestep_count)
379                   ENDDO
380                 
381                   sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
382 
383                ENDDO
384            ENDDO
385
386          CASE ( 'v' )
387
388             DO  i = nxl, nxr
389                DO  j = nys, nyn
390
391                   DO  k = nzb_u_inner(j,i)+1, nzt
392
393                      tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +     &
394                                     vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
395
396                      tend(k,j,i) = tend(k,j,i) + tmp_tend
397
398                      sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend *             &
399                                     weight_substep(intermediate_timestep_count)
400                   ENDDO
401                 
402                   sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
403
404                ENDDO
405            ENDDO
406
407          CASE ( 'pt' )
408
409             DO  i = nxl, nxr
410                DO  j = nys, nyn
411
412                   DO  k = nzb_u_inner(j,i)+1, nzt
413
414                      tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +    &
415                                     ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
416
417                      tend(k,j,i) = tend(k,j,i) + tmp_tend
418
419                      sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend *             &
420                                     weight_substep(intermediate_timestep_count)
421                   ENDDO
422
423                   sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
424
425                ENDDO
426            ENDDO
427
428          CASE ( 'q' )
429
430             DO  i = nxl, nxr
431                DO  j = nys, nyn
432
433                   DO  k = nzb_u_inner(j,i)+1, nzt
434
435                      tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +    &
436                                     qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
437
438                      tend(k,j,i) = tend(k,j,i) + tmp_tend
439
440                      sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend *             &
441                                     weight_substep(intermediate_timestep_count)
442                   ENDDO
443                 
444                   sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
445
446                ENDDO
447            ENDDO
448
449          CASE DEFAULT
450             message_string = 'unknown prognostic variable "' // prog_var // '"'
451             CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
452
453       END SELECT
454
455    END SUBROUTINE nudge
456
457
458!------------------------------------------------------------------------------!
459! Description:
460! ------------
461!> Call for grid point i,j
462!------------------------------------------------------------------------------!
463
464    SUBROUTINE nudge_ij( i, j, time, prog_var )
465
466       USE arrays_3d,                                                          &
467           ONLY:  pt, ptnudge, q, qnudge, tend, timenudge, tmp_tnudge, tnudge, &
468                  u, unudge, v, vnudge
469
470       USE control_parameters,                                                 &
471           ONLY:  dt_3d, intermediate_timestep_count, message_string
472
473       USE indices,                                                            &
474           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_u_inner, nzt
475
476       USE kinds
477
478       USE statistics,                                                         &
479           ONLY:  hom, sums_ls_l, weight_substep
480
481       IMPLICIT NONE
482
483
484       CHARACTER (LEN=*) ::  prog_var  !<
485
486       REAL(wp) ::  tmp_tend    !<
487       REAL(wp) ::  dtm         !<
488       REAL(wp) ::  dtp         !<
489       REAL(wp) ::  time        !<
490
491       INTEGER(iwp) ::  i  !<
492       INTEGER(iwp) ::  j  !<
493       INTEGER(iwp) ::  k  !<
494       INTEGER(iwp) ::  nt  !<
495
496
497       nt = 1
498       DO WHILE ( time > timenudge(nt) )
499         nt = nt+1
500       ENDDO
501       IF ( time /= timenudge(1) )  THEN
502         nt = nt-1
503       ENDIF
504
505       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
506       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
507
508       SELECT CASE ( prog_var )
509
510          CASE ( 'u' )
511
512             DO  k = nzb_u_inner(j,i)+1, nzt
513
514                tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +           &
515                               unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
516
517                tend(k,j,i) = tend(k,j,i) + tmp_tend
518
519                sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend                     &
520                                 * weight_substep(intermediate_timestep_count)
521             ENDDO
522
523             sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
524
525          CASE ( 'v' )
526
527             DO  k = nzb_u_inner(j,i)+1, nzt
528
529                tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +           &
530                               vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
531
532                tend(k,j,i) = tend(k,j,i) + tmp_tend
533
534                sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend                     &
535                                 * weight_substep(intermediate_timestep_count)
536             ENDDO
537
538             sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
539
540          CASE ( 'pt' )
541
542             DO  k = nzb_u_inner(j,i)+1, nzt
543
544                tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +          &
545                               ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
546
547                tend(k,j,i) = tend(k,j,i) + tmp_tend
548
549                sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend                     &
550                                 * weight_substep(intermediate_timestep_count)
551             ENDDO
552
553             sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
554
555
556          CASE ( 'q' )
557
558             DO  k = nzb_u_inner(j,i)+1, nzt
559
560                tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +          &
561                               qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
562
563                tend(k,j,i) = tend(k,j,i) + tmp_tend
564
565                sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend                     &
566                                 * weight_substep(intermediate_timestep_count)
567             ENDDO
568
569             sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
570
571          CASE DEFAULT
572             message_string = 'unknown prognostic variable "' // prog_var // '"'
573             CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
574
575       END SELECT
576
577
578    END SUBROUTINE nudge_ij
579
580
581!------------------------------------------------------------------------------!
582! Description:
583! ------------
584!> @todo Missing subroutine description.
585!------------------------------------------------------------------------------!
586    SUBROUTINE nudge_ref ( time )
587
588       USE arrays_3d,                                                          &
589           ONLY:  time_vert, ptnudge, pt_init, qnudge, q_init, unudge, u_init, &
590                  vnudge, v_init
591
592       USE kinds
593
594
595       IMPLICIT NONE
596
597       INTEGER(iwp) ::  nt                    !<
598
599       REAL(wp)             ::  fac           !<
600       REAL(wp), INTENT(in) ::  time          !<
601
602!
603!--    Interpolation in time of NUDGING_DATA for pt_init and q_init. This is
604!--    needed for correct upper boundary conditions for pt and q and in case that
605!      large scale subsidence as well as scalar Rayleigh-damping are used
606       nt = 1
607       DO WHILE ( time > time_vert(nt) )
608          nt = nt + 1
609       ENDDO
610       IF ( time /= time_vert(nt) )  THEN
611        nt = nt - 1
612       ENDIF
613
614       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
615
616       pt_init = ptnudge(:,nt) + fac * ( ptnudge(:,nt+1) - ptnudge(:,nt) )
617       q_init  = qnudge(:,nt) + fac * ( qnudge(:,nt+1) - qnudge(:,nt) )
618       u_init  = unudge(:,nt) + fac * ( unudge(:,nt+1) - unudge(:,nt) )
619       v_init  = vnudge(:,nt) + fac * ( vnudge(:,nt+1) - vnudge(:,nt) )
620
621    END SUBROUTINE nudge_ref
622
623 END MODULE nudge_mod
Note: See TracBrowser for help on using the repository browser.