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

Last change on this file since 1850 was 1850, checked in by maronga, 6 years ago

added _mod string to several filenames to meet the naming convection for modules

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