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

Last change on this file since 1321 was 1321, checked in by raasch, 11 years ago

last commit documented

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