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

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

adjustments in case of nudging + bugfix for KIND in CMPLX function

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