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

Last change on this file since 1682 was 1682, checked in by knoop, 8 years ago

Code annotations made doxygen readable

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