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

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

last commit documented

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