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

Last change on this file since 1380 was 1380, checked in by heinze, 10 years ago

Upper boundary conditions for pt and q in case of nudging adjusted

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