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

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

large scale forcing enabled

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