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

Last change on this file since 1366 was 1366, checked in by boeske, 8 years ago

last commit documented

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