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

Last change on this file since 1384 was 1383, checked in by boeske, 10 years ago

last commit documented

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