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

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

last commit documented

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