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

Last change on this file since 1242 was 1242, checked in by heinze, 8 years ago

Last commmit documented

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