source: palm/trunk/SOURCE/nudging.f90 @ 1748

Last change on this file since 1748 was 1683, checked in by knoop, 8 years ago

last commit documented

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