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

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

last commit documented / copyright update

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