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

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

Error message specified.

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