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

Last change on this file since 1413 was 1399, checked in by heinze, 11 years ago

last commit documented

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