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

Last change on this file since 1357 was 1356, checked in by heinze, 11 years ago

last commit documented

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