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

Last change on this file since 1757 was 1757, checked in by maronga, 6 years ago

some changes in land surface model, radiation model, nudging and some minor updates

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