Newer
Older
WAUTELET Philippe
committed
!MNH_LIC Copyright 2006-2023 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
!-----------------------------------------------------------------
! ####################
MODULE MODI_ZDIFFUSETUP
! ####################
!
INTERFACE
!
! ######################################################################
SUBROUTINE ZDIFFUSETUP (PZZ,&
PZDIFFU_HALO2)
! ######################################################################
!JUAN
USE MODE_TYPE_ZDIFFU
IMPLICIT NONE
!JUAN
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height
!JUAN
TYPE(TYPE_ZDIFFU_HALO2) :: PZDIFFU_HALO2
!JUAN
END SUBROUTINE ZDIFFUSETUP
!
END INTERFACE
!
END MODULE MODI_ZDIFFUSETUP
!
!
!
! #########################################################################
SUBROUTINE ZDIFFUSETUP (PZZ,&
PZDIFFU_HALO2)
! #########################################################################
!
!!**** *ZDIFFUSETUP* - routine to calculate the interpolation coefficient needed
!! for the truly horizontal diffusion scheme
!!
!! REFERENCE
!! ---------
!!
WAUTELET Philippe
committed
!! Zangl, G., 2002: An improved method for computing horizontal diffusion in a
!! sigma-coordinate model and its application to simulations
!! over mountainous topography. Mon. Wea. Rev. 130, 1423-1432.
!!
!! AUTHOR
!! ------
!!
WAUTELET Philippe
committed
!! G. Zangl * University of Munich*
!
! Modifications:
! J. Escobar 07/10/2015: remove print
! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
WAUTELET Philippe
committed
! P. Wautelet 26/04/2019: use modd_precision parameters for datatypes of MPI communications
WAUTELET Philippe
committed
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
!
!* 0. DECLARATIONS
! ------------
WAUTELET Philippe
committed
USE MODD_ARGSLIST_ll, ONLY: LIST_ll, HALO2LIST_ll
USE MODD_CONF
WAUTELET Philippe
committed
use modd_precision, only: MNHINT_MPI
USE MODD_VAR_ll, ONLY: NMNH_COMM_WORLD
!
WAUTELET Philippe
committed
USE MODE_SUM_LL
USE MODE_TYPE_ZDIFFU
WAUTELET Philippe
committed
USE MODI_RELAX
USE MODI_SHUMAN
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
!
IMPLICIT NONE
!
!* 0.1 declarations of arguments
!
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height
!JUAN
TYPE(TYPE_ZDIFFU_HALO2) :: PZDIFFU_HALO2
!JUAN
! local variables
INTEGER :: IIB_ll,IJB_ll ! Lower bounds of the physical
! global domain in x and y directions
INTEGER :: IIE_ll,IJE_ll ! Upper bounds of the physical
! gloBal domain in x and y directions
! sub-domain:
INTEGER :: IIB,IJB,IKB ! Lower bounds of the physical
! sub-domain in x,y and z directions
INTEGER :: IIE,IJE,IKE ! Upper bounds of the physical
! sub-domain in x,y and z directions
INTEGER :: IIU,IJU,IKU ! domain sizes
INTEGER :: JI,JJ,JK ! Loop indexes
INTEGER :: IIMAX_ll,IJMAX_ll ! Number of points of
! Global physical domain
! in the x and y directions
INTEGER, DIMENSION(:,:), ALLOCATABLE :: IKIP1,IKIP2,IKIM1,IKIM2,IKJP1,IKJP2,IKJM1,IKJM2,IKMAX
REAL, DIMENSION(:,:), ALLOCATABLE :: ZN4HGTI,ZN4HGTJ,ZMDHGTI,ZMDHGTJ
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZZMASS
REAL SCAL
!JUAN
INTEGER, DIMENSION(:,:), ALLOCATABLE :: IKIP1_HALO2,IKIP2_HALO2,IKIM1_HALO2,IKIM2_HALO2
INTEGER, DIMENSION(:,:), ALLOCATABLE :: IKJP1_HALO2,IKJP2_HALO2,IKJM1_HALO2,IKJM2_HALO2
INTEGER, DIMENSION(:,:), ALLOCATABLE :: IKMAX_HALO2
REAL, DIMENSION(:,:), ALLOCATABLE :: ZN4HGTI_HALO2,ZN4HGTJ_HALO2,ZMDHGTI_HALO2,ZMDHGTJ_HALO2
TYPE(LIST_ll), POINTER :: TZHGTMASS_ll
TYPE(HALO2LIST_ll), POINTER :: TZHGTHALO2_ll
INTEGER :: KZDLB_ll,IERR
!JUAN
INTEGER:: IERROR ! DUMMY VARIABLE FOR ERROR MESSAGES
!
!* 1. COMPUTE THE PHYSICAL SUBDOMAIN BOUNDS
! ---------------------------------------
!
CALL GET_INDICE_ll( IIB,IJB,IIE,IJE)
IKE = SIZE(PZZ,3) - JPVEXT
IKB = 1 + JPVEXT
! Global physical dimensions
CALL GET_GLOBALDIMS_ll ( IIMAX_ll,IJMAX_ll)
IIU = SIZE(PZZ,1)
IJU = SIZE(PZZ,2)
IKU = SIZE(PZZ,3)
!PRINT*,'Interpolation coefficients for truly horizontal diffusion are computed'
ALLOCATE (ZZMASS(IIU,IJU,IKU))
!JUAN
ALLOCATE (IKIP1_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),IKIP2_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),&
IKIM1_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),IKIM2_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),&
IKJP1_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),IKJP2_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),&
IKJM1_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),IKJM2_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),&
IKMAX_HALO2(IIB-2:IIE+2,IJB-2:IJE+2))
ALLOCATE (ZN4HGTI_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),ZN4HGTJ_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),&
ZMDHGTI_HALO2(IIB-2:IIE+2,IJB-2:IJE+2),ZMDHGTJ_HALO2(IIB-2:IIE+2,IJB-2:IJE+2))
!JUAN
NULLIFY(TZHGTMASS_ll,TZHGTHALO2_ll)
! Compute height field at mass points
WAUTELET Philippe
committed
ZZMASS = MZF(PZZ)
WAUTELET Philippe
committed
CALL INIT_HALO2_ll( TZHGTHALO2_ll, 1, IIU, IJU, IKU, 'ZDIFFSETUP::HGTHALO2' )
WAUTELET Philippe
committed
CALL ADD3DFIELD_ll( TZHGTMASS_ll, ZZMASS, 'ZDIFFUSETUP::ZZMASS' )
CALL UPDATE_HALO2_ll(TZHGTMASS_ll,TZHGTHALO2_ll,IERROR)
!JUAN
PZDIFFU_HALO2%XZZ(IIB-2,IJB:IJE,:) = TZHGTHALO2_ll%HALO2%WEST(IJB:IJE,:)
PZDIFFU_HALO2%XZZ(IIE+2,IJB:IJE,:) = TZHGTHALO2_ll%HALO2%EAST(IJB:IJE,:)
PZDIFFU_HALO2%XZZ(IIB:IIE,IJB-2,:) = TZHGTHALO2_ll%HALO2%SOUTH(IIB:IIE,:)
PZDIFFU_HALO2%XZZ(IIB:IIE,IJE+2,:) = TZHGTHALO2_ll%HALO2%NORTH(IIB:IIE,:)
PZDIFFU_HALO2%XZZ(1:IIU,1:IJU,1:IKU) = ZZMASS
!print *,"zdiffu :: IIB-2",IIB-2
!print *,"zdiffu :: IIE+2",IIE+2
!print *,"zdiffu :: IJB-2",IJB-2
!print *,"zdiffu :: IJE+2",IJE+2
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
!DO JI=IIB-2,IIE+2
!print *,"zdiffu :: PZDIFFU_HALO2%X=JI=",JI,(PZDIFFU_HALO2%XZZ(JI,JJ,2),JJ=IJB-2,IJE+2)
!ENDDO
! Compute the vertical index of the remote grid points having the same height as the local
! grid point. This index is a real number composed of the index of the lower neighbouring model
! level and the weighting coefficient needed for the linear vertical interpolation between the
! model levels.
!JUAN
!JUAN
!JUAN
IKIP1_HALO2 = IKB
IKIP2_HALO2 = IKB
IKIM1_HALO2 = IKB
IKIM2_HALO2 = IKB
IKJP1_HALO2 = IKB
IKJP2_HALO2 = IKB
IKJM1_HALO2 = IKB
IKJM2_HALO2 = IKB
PZDIFFU_HALO2%XRKIP1 = IKB
PZDIFFU_HALO2%XRKIP2 = IKB
PZDIFFU_HALO2%XRKIM1 = IKB
PZDIFFU_HALO2%XRKIM2 = IKB
PZDIFFU_HALO2%XRKJP1 = IKB
PZDIFFU_HALO2%XRKJP2 = IKB
PZDIFFU_HALO2%XRKJM1 = IKB
PZDIFFU_HALO2%XRKJM2 = IKB
PZDIFFU_HALO2%XREDFACI = 1.0
PZDIFFU_HALO2%XREDFACJ = 1.0
CALL INDINT_HALO2( 1, 0,PZDIFFU_HALO2%XZZ,PZDIFFU_HALO2%XRKIP1,IKIP1_HALO2,IIB,IJB)
CALL INDINT_HALO2( 2, 0,PZDIFFU_HALO2%XZZ,PZDIFFU_HALO2%XRKIP2,IKIP2_HALO2,IIB,IJB)
CALL INDINT_HALO2(-1, 0,PZDIFFU_HALO2%XZZ,PZDIFFU_HALO2%XRKIM1,IKIM1_HALO2,IIB,IJB)
CALL INDINT_HALO2(-2, 0,PZDIFFU_HALO2%XZZ,PZDIFFU_HALO2%XRKIM2,IKIM2_HALO2,IIB,IJB)
CALL INDINT_HALO2( 0, 1,PZDIFFU_HALO2%XZZ,PZDIFFU_HALO2%XRKJP1,IKJP1_HALO2,IIB,IJB)
CALL INDINT_HALO2( 0, 2,PZDIFFU_HALO2%XZZ,PZDIFFU_HALO2%XRKJP2,IKJP2_HALO2,IIB,IJB)
CALL INDINT_HALO2( 0,-1,PZDIFFU_HALO2%XZZ,PZDIFFU_HALO2%XRKJM1,IKJM1_HALO2,IIB,IJB)
CALL INDINT_HALO2( 0,-2,PZDIFFU_HALO2%XZZ,PZDIFFU_HALO2%XRKJM2,IKJM2_HALO2,IIB,IJB)
PZDIFFU_HALO2%NZDI = MAX(IKIP1_HALO2,IKIP2_HALO2,IKIM1_HALO2,IKIM2_HALO2)
PZDIFFU_HALO2%NZDJ = MAX(IKJP1_HALO2,IKJP2_HALO2,IKJM1_HALO2,IKJM2_HALO2)
IKMAX_HALO2 = MAX(PZDIFFU_HALO2%NZDI,PZDIFFU_HALO2%NZDJ)
PZDIFFU_HALO2%NZDLB = MAXVAL(IKMAX_HALO2) ! Model level, above which a truly horizontal computation of diffusion
! is possible at all grid points
!JUAN
WAUTELET Philippe
committed
CALL MPI_ALLREDUCE(PZDIFFU_HALO2%NZDLB ,KZDLB_ll, 1, MNHINT_MPI, MPI_MAX, NMNH_COMM_WORLD, IERR)
!print*,"zdiffusetup:: PZDIFFU_HALO2%NZDLB=",PZDIFFU_HALO2%NZDLB,KZDLB_ll
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
PZDIFFU_HALO2%NZDLB = KZDLB_ll
!JUAN
! Compute reduction factors for diffusion coefficient in I and J direction
! Their purpose is to decrease the diffusion coefficient when there is a large height difference
! among the grid points entering into the computation of diffusion
SCAL = 1.E9 ! scaling factor for coefficient reduction
DO JI = IIB,IIE
DO JJ = IJB-1,IJE+1
ZN4HGTI_HALO2(JI,JJ) = PZDIFFU_HALO2%XZZ(JI+2,JJ,IKB)+PZDIFFU_HALO2%XZZ(JI-2,JJ,IKB)- &
4*(PZDIFFU_HALO2%XZZ(JI+1,JJ,IKB)+PZDIFFU_HALO2%XZZ(JI-1,JJ,IKB))+ &
6* PZDIFFU_HALO2%XZZ(JI,JJ,IKB)
ZMDHGTI_HALO2(JI,JJ) = (PZDIFFU_HALO2%XZZ(JI+2,JJ,IKB)+PZDIFFU_HALO2%XZZ(JI-2,JJ,IKB)+ &
PZDIFFU_HALO2%XZZ(JI+1,JJ,IKB)+PZDIFFU_HALO2%XZZ(JI-1,JJ,IKB)- &
4* PZDIFFU_HALO2%XZZ(JI,JJ,IKB))/4.
PZDIFFU_HALO2%XREDFACI(JI,JJ) = SCAL/(SCAL+ZN4HGTI_HALO2(JI,JJ)**4+ZMDHGTI_HALO2(JI,JJ)**4)
IF (((ZN4HGTI_HALO2(JI,JJ).GT.100).AND.(ZMDHGTI_HALO2(JI,JJ).GT.10)).OR. &
((ZN4HGTI_HALO2(JI,JJ).GT.10).AND.(ZMDHGTI_HALO2(JI,JJ).GT.100))) THEN
PZDIFFU_HALO2%XREDFACI(JI,JJ) = MIN(PZDIFFU_HALO2%XREDFACI(JI,JJ),0.1)
ENDIF
ENDDO
ENDDO
DO JI = IIB-1,IIE+1
DO JJ = IJB,IJE
ZN4HGTJ_HALO2(JI,JJ) = PZDIFFU_HALO2%XZZ(JI,JJ+2,IKB)+PZDIFFU_HALO2%XZZ(JI,JJ-2,IKB)- &
4*(PZDIFFU_HALO2%XZZ(JI,JJ+1,IKB)+PZDIFFU_HALO2%XZZ(JI,JJ-1,IKB))+ &
6* PZDIFFU_HALO2%XZZ(JI,JJ,IKB)
ZMDHGTJ_HALO2(JI,JJ) = (PZDIFFU_HALO2%XZZ(JI,JJ+2,IKB)+PZDIFFU_HALO2%XZZ(JI,JJ-2,IKB)+ &
PZDIFFU_HALO2%XZZ(JI,JJ+1,IKB)+PZDIFFU_HALO2%XZZ(JI,JJ-1,IKB)- &
4* PZDIFFU_HALO2%XZZ(JI,JJ,IKB))/4.
PZDIFFU_HALO2%XREDFACJ(JI,JJ) = SCAL/(SCAL+ZN4HGTJ_HALO2(JI,JJ)**4+ZMDHGTJ_HALO2(JI,JJ)**4)
IF (((ZN4HGTJ_HALO2(JI,JJ).GT.100).AND.(ZMDHGTJ_HALO2(JI,JJ).GT.10)).OR. &
((ZN4HGTJ_HALO2(JI,JJ).GT.10).AND.(ZMDHGTJ_HALO2(JI,JJ).GT.100))) THEN
PZDIFFU_HALO2%XREDFACJ(JI,JJ) = MIN(PZDIFFU_HALO2%XREDFACJ(JI,JJ),0.1)
ENDIF
ENDDO
ENDDO
CALL DEL_HALO2_ll( TZHGTHALO2_ll )
CONTAINS
SUBROUTINE INDINT_HALO2(KII,KIJ,PZMASS,PKIND,KKMIN,KIB,KJB)
use mode_msg
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
IMPLICIT NONE
INTEGER, INTENT(IN) :: KII,KIJ ! Relative position of remote points
INTEGER, INTENT(IN) :: KIB,KJB ! definition of domain begin
REAL, DIMENSION(KIB-2:,KJB-2:,:), INTENT(IN) :: PZMASS ! Height of mass points
REAL, DIMENSION(KIB-2:,KJB-2:,:), INTENT(INOUT) :: PKIND ! Real k index for vertical interpolation
INTEGER, DIMENSION(KIB-2:,KJB-2:),INTENT(INOUT) :: KKMIN ! Lowest model level for which truly horizontal computation of
! the diffusion is possible without intersecting the ground
! Local variables
! Domain sizes
INTEGER IIB,IJB,IIE,IJE,IKB,IKE,II1,II2,IJ1,IJ2
! Loop indices
INTEGER JI,JJ,JK,JIR,JJR,JK2
REAL ZHGT ! Height of the local grid point
CALL GET_INDICE_ll( IIB,IJB,IIE,IJE)
IKE = SIZE(PZMASS,3) - JPVEXT
IKB = 1 + JPVEXT ! (is this correct?)
IF ((KII.EQ.0).AND.(KIJ.NE.0)) THEN
II1 = IIB-1
II2 = IIE+1
!JUAN II1 = IIB-2
!JUAN II2 = IIE+2
IJ1 = IJB
IJ2 = IJE
ELSE IF ((KIJ.EQ.0).AND.(KII.NE.0)) THEN
II1 = IIB
II2 = IIE
IJ1 = IJB-1
IJ2 = IJE+1
!JUAN IJ1 = IJB-2
!JUAN IJ2 = IJE+2
ELSE
call Print_msg( NVERB_FATAL, 'GEN', 'INDINT_HALO2', 'KII=0 and KIJ=0' )
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
ENDIF
DO JI = II1,II2
DO JJ = IJ1,IJ2
KKMIN(JI,JJ) = IKB
DO JK = IKB,IKE
JIR = JI + KII ! location of the remote points for which the real index is to be computed
JJR = JJ + KIJ
ZHGT = PZMASS(JI,JJ,JK) ! height for which the index is needed
IF (ZHGT.GT.PZMASS(JIR,JJR,JK)) THEN ! local point higher than remote point; search upward
IF (ZHGT.GT.PZMASS(JIR,JJR,IKE)) THEN
PKIND(JI,JJ,JK) = IKE
GOTO 20
ENDIF
DO JK2 = JK,IKE-1
IF((ZHGT.GT.PZMASS(JIR,JJR,JK2)).AND.(ZHGT.LE.PZMASS(JIR,JJR,JK2+1))) THEN
PKIND(JI,JJ,JK) = JK2+(ZHGT-PZMASS(JIR,JJR,JK2))/&
(PZMASS(JIR,JJR,JK2+1)-PZMASS(JIR,JJR,JK2))
GOTO 20
ENDIF
ENDDO
20 CONTINUE
ELSE ! local point lower than remote point; search downward
IF (ZHGT.LT.PZMASS(JIR,JJR,IKB)) THEN
PKIND(JI,JJ,JK) = IKB
KKMIN(JI,JJ) = MAX(KKMIN(JI,JJ),JK+1)
GOTO 25
ELSE IF (ZHGT.EQ.PZMASS(JIR,JJR,IKB)) THEN
PKIND(JI,JJ,JK) = IKB
KKMIN(JI,JJ) = MAX(KKMIN(JI,JJ),JK)
GOTO 25
ENDIF
DO JK2 = JK,IKB+1,-1
IF((ZHGT.GT.PZMASS(JIR,JJR,JK2-1)).AND.(ZHGT.LE.PZMASS(JIR,JJR,JK2))) THEN
PKIND(JI,JJ,JK) = JK2-1+(ZHGT-PZMASS(JIR,JJR,JK2-1))/&
(PZMASS(JIR,JJR,JK2)-PZMASS(JIR,JJR,JK2-1))
GOTO 25
ENDIF
ENDDO
25 CONTINUE
ENDIF
ENDDO
ENDDO
ENDDO
IF ( MINVAL(KKMIN) == 0 ) THEN
call Print_msg( NVERB_FATAL, 'GEN', 'INDINT_HALO2', 'MINVAL(KKMIN)=0' )
IF ( MINVAL(INT(PKIND)) == 0 ) THEN
call Print_msg( NVERB_FATAL, 'GEN', 'INDINT_HALO2', 'MINVAL(INT(PKIND))=0' )
END SUBROUTINE INDINT_HALO2
END SUBROUTINE ZDIFFUSETUP