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

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

minor changes in profile data output of lsf tendencies, variables renamed

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