source: palm/trunk/SOURCE/nudging_mod.f90 @ 2000

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

Forced header and separation lines into 80 columns

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