Changeset 690 for trunk/NEMO/OPA_SRC
- Timestamp:
- 2007-07-04T15:24:29+02:00 (17 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90
r689 r690 1 1 !!---------------------------------------------------------------------- 2 !! *** ldfdyn_c2d.h90 *** 3 !!---------------------------------------------------------------------- 4 !! ldf_dyn_c2d : set the lateral viscosity coefficients 5 !! ldf_dyn_c2d_orca : specific case for orca r2 and r4 6 !!---------------------------------------------------------------------- 7 8 !!---------------------------------------------------------------------- 9 !! OPA 9.0 , LOCEAN-IPSL (2005) 10 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90,v 1.9 2007/06/29 17:01:51 opalod Exp $ 2 !! *** ldfdyn_c3d.h90 *** 3 !!---------------------------------------------------------------------- 4 5 !!---------------------------------------------------------------------- 6 !! OPA 9.0 , LOCEAN-IPSL (2005) 7 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90,v 1.9 2007/07/04 13:02:29 opalod Exp $ 11 8 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 12 9 !!---------------------------------------------------------------------- 13 10 14 SUBROUTINE ldf_dyn_c2d( ld_print ) 15 !!---------------------------------------------------------------------- 16 !! *** ROUTINE ldf_dyn_c2d *** 17 !! 11 !!---------------------------------------------------------------------- 12 !! 'key_dynldf_c3d' 3D lateral eddy viscosity coefficients 13 !!---------------------------------------------------------------------- 14 15 SUBROUTINE ldf_dyn_c3d( ld_print ) 16 !!---------------------------------------------------------------------- 17 !! *** ROUTINE ldf_dyn_c3d *** 18 !! 18 19 !! ** Purpose : initializations of the horizontal ocean physics 19 20 !! 20 !! ** Method : 21 !! 2D eddy viscosity coefficients ( longitude, latitude ) 22 !! 23 !! harmonic operator : ahm1 is defined at t-point 24 !! ahm2 is defined at f-point 25 !! + isopycnal : ahm3 is defined at u-point 26 !! or geopotential ahm4 is defined at v-point 27 !! iso-model level : ahm3, ahm4 not used 28 !! 29 !! biharmonic operator : ahm1 is defined at u-point 30 !! ahm2 is defined at v-point 31 !! : ahm3, ahm4 not used 32 !! 33 !!---------------------------------------------------------------------- 21 !! ** Method : 3D eddy viscosity coef. ( longitude, latitude, depth ) 22 !! laplacian operator : ahm1, ahm2 defined at T- and F-points 23 !! ahm2, ahm4 never used 24 !! bilaplacian operator : ahm1, ahm2 never used 25 !! : ahm3, ahm4 defined at U- and V-points 26 !! ??? explanation of the default is missing 27 !!---------------------------------------------------------------------- 28 !! * Modules used 29 USE ldftra_oce, ONLY : aht0 30 34 31 !! * Arguments 35 32 LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout 36 33 37 !! * Local variables 38 REAL(wp) :: za00, zdx_max 34 !! * local variables 35 INTEGER :: ji, jj, jk ! dummy loop indices 36 REAL(wp) :: & 37 zr = 0.2 , & ! maximum of the reduction factor at the bottom ocean 38 ! ! ( 0 < zr < 1 ) 39 zh = 500., & ! depth of at which start the reduction ( > dept(1) ) 40 zdx_max , & ! maximum grid spacing over the global domain 41 za00, zc, zd ! temporary scalars 42 REAL(wp), DIMENSION(jpk) :: zcoef ! temporary workspace 39 43 !!---------------------------------------------------------------------- 40 44 41 45 IF(lwp) WRITE(numout,*) 42 IF(lwp) WRITE(numout,*) 'ldf_dyn_c 2d : 2dlateral eddy viscosity coefficient'46 IF(lwp) WRITE(numout,*) 'ldf_dyn_c3d : 3D lateral eddy viscosity coefficient' 43 47 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 44 IF(lwp) WRITE(numout,*) 45 46 ! harmonic operator (ahm1, ahm2) :( T- and F- points) (used for laplacian operators47 ! ================= ==============whatever its orientation is)48 49 50 ! Set ahm1 and ahm2 ( T- and F- points) (used for laplacian operators 51 ! ================= whatever its orientation is) 48 52 IF( ln_dynldf_lap ) THEN 49 53 ! define ahm1 and ahm2 at the right grid point position … … 57 61 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0 58 62 63 59 64 za00 = ahm0 / zdx_max 60 ahm1(:,:) = za00 * e1t(:,:)61 ahm2(:,:) = za00 * e1f(:,:)62 65 63 66 IF( ln_dynldf_iso ) THEN … … 68 71 ENDIF 69 72 73 CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm1 ) ! vertical profile 74 CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm2 ) ! vertical profile 75 DO jk = 1,jpk 76 ahm1(:,:,jk) = za00 * e1t(:,:) * ahm1(:,:,jk) 77 ahm2(:,:,jk) = za00 * e1f(:,:) * ahm2(:,:,jk) 78 END DO 79 80 70 81 ! Special case for ORCA R2 and R4 configurations (overwrite the value of ahm1 ahm2) 71 82 ! ============================================== 72 IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) ) CALL ldf_dyn_c2d_orca( ld_print ) 83 IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) ) THEN 84 IF(lwp) WRITE(numout,*) 85 IF(lwp) WRITE(numout,*) ' ORCA R2 or R4: overwrite the previous definition of ahm' 86 IF(lwp) WRITE(numout,*) ' =============' 87 CALL ldf_dyn_c3d_orca( ld_print ) 88 ENDIF 89 90 ENDIF 91 92 ! Control print 93 IF(lwp .AND. ld_print ) THEN 94 WRITE(numout,*) 95 WRITE(numout,*) ' 3D ahm1 array (k=1)' 96 CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout ) 97 WRITE(numout,*) 98 WRITE(numout,*) ' 3D ahm2 array (k=1)' 99 CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout ) 100 ENDIF 101 102 103 ! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator 104 ! ================================ whatever its orientation is) 105 ! (USER: modify ahm3 and ahm4 following your desiderata) 106 ! Here: ahm is proportional to the cube of the maximum of the gridspacing 107 ! in the to horizontal direction 108 109 IF( ln_dynldf_bilap ) THEN 110 111 zdx_max = MAXVAL( e1u(:,:) ) 112 IF( lk_mpp ) CALL mpp_max( zdx_max ) ! max over the global domain 113 114 IF(lwp) WRITE(numout,*) ' bi-laplacian operator: ahm proportional to e1**3 ' 115 IF(lwp) WRITE(numout,*) ' Caution, here we assume your mesh is isotropic ...' 116 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0 117 118 za00 = ahm0 / ( zdx_max * zdx_max * zdx_max ) 119 ahm3(:,:,1) = za00 * e1u(:,:) * e1u(:,:) * e1u(:,:) 120 ahm4(:,:,1) = za00 * e1v(:,:) * e1v(:,:) * e1v(:,:) 121 122 zh = MAX( zh, fsdept(1,1,1) ) ! at least the first reach ahm0 123 IF( ln_zco ) THEN ! z-coordinate, same profile everywhere 124 IF(lwp) WRITE(numout,'(36x," ahm ", 7x)') 125 DO jk = 1, jpk 126 IF( fsdept(1,1,jk) <= zh ) THEN 127 zcoef(jk) = 1.e0 128 ELSE 129 zcoef(jk) = 1.e0 + ( zr - 1.e0 ) & 130 & * ( 1. - EXP( ( fsdept(1,1,jk ) - zh ) / zh ) ) & 131 & / ( 1. - EXP( ( fsdept(1,1,jpkm1) - zh ) / zh ) ) 132 ENDIF 133 ahm3(:,:,jk) = ahm3(:,:,1) * zcoef(jk) 134 ahm4(:,:,jk) = ahm4(:,:,1) * zcoef(jk) 135 IF(lwp) WRITE(numout,'(34x,E7.2,8x,i3)') zcoef(jk) * ahm0, jk 136 END DO 137 ELSE ! partial steps or s-ccordinate 138 # if defined key_zco 139 zc = gdept_0(jpkm1) 140 # else 141 zc = MAXVAL( fsdept(:,:,jpkm1) ) 142 # endif 143 IF( lk_mpp ) CALL mpp_max( zc ) ! max over the global domain 144 145 zc = 1. / ( 1. - EXP( ( zc - zh ) / zh ) ) 146 DO jk = 2, jpkm1 147 DO jj = 1, jpj 148 DO ji = 1, jpi 149 IF( fsdept(ji,jj,jk) <= zh ) THEN 150 ahm3(ji,jj,jk) = ahm3(ji,jj,1) 151 ahm4(ji,jj,jk) = ahm4(ji,jj,1) 152 ELSE 153 zd = 1.e0 + ( zr - 1.e0 ) * ( 1. - EXP( ( fsdept(ji,jj,jk) - zh ) / zh ) ) * zc 154 ahm3(ji,jj,jk) = ahm3(ji,jj,1) * zd 155 ahm4(ji,jj,jk) = ahm4(ji,jj,1) * zd 156 ENDIF 157 END DO 158 END DO 159 END DO 160 ahm3(:,:,jpk) = ahm3(:,:,jpkm1) 161 ahm4(:,:,jpk) = ahm4(:,:,jpkm1) 162 IF(lwp) WRITE(numout,'(36x," ahm ", 7x)') 163 DO jk = 1, jpk 164 IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(1,1,jk), jk 165 END DO 166 ENDIF 73 167 74 168 ! Control print 75 169 IF( lwp .AND. ld_print ) THEN 76 170 WRITE(numout,*) 77 WRITE(numout,*) 'inildf: 2D ahm1 array'78 CALL prihre(ahm 1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)171 WRITE(numout,*) 'inildf: ahm3 array at level 1' 172 CALL prihre(ahm3(:,:,1 ),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 79 173 WRITE(numout,*) 80 WRITE(numout,*) 'inildf: 2D ahm2 array'81 CALL prihre(ahm 2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)174 WRITE(numout,*) 'inildf: ahm4 array at level 1' 175 CALL prihre(ahm4(:,:,jpk),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 82 176 ENDIF 83 177 ENDIF 84 178 85 ! biharmonic operator (ahm3, ahm4) : at U- and V-points (used for bilaplacian operator 86 ! ================================= whatever its orientation is) 87 IF( ln_dynldf_bilap ) THEN 88 ! (USER: modify ahm3 and ahm4 following your desiderata) 89 ! Here: ahm is proportional to the cube of the maximum of the gridspacing 90 ! in the to horizontal direction 91 92 zdx_max = MAXVAL( e1u(:,:) ) 93 IF( lk_mpp ) CALL mpp_max( zdx_max ) ! max over the global domain 94 95 IF(lwp) WRITE(numout,*) ' bi-laplacian operator: ahm proportional to e1**3 ' 96 IF(lwp) WRITE(numout,*) ' Caution, here we assume your mesh is isotropic ...' 97 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0 98 99 za00 = ahm0 / ( zdx_max * zdx_max * zdx_max ) 100 ahm3(:,:) = za00 * e1u(:,:) * e1u(:,:) * e1u(:,:) 101 ahm4(:,:) = za00 * e1v(:,:) * e1v(:,:) * e1v(:,:) 102 103 ! Control print 104 IF( lwp .AND. ld_print ) THEN 105 WRITE(numout,*) 106 WRITE(numout,*) 'inildf: ahm3 array' 107 CALL prihre(ahm3,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 108 WRITE(numout,*) 109 WRITE(numout,*) 'inildf: ahm4 array' 110 CALL prihre(ahm4,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 111 ENDIF 112 ENDIF 113 114 115 END SUBROUTINE ldf_dyn_c2d 116 117 118 SUBROUTINE ldf_dyn_c2d_orca( ld_print ) 119 !!---------------------------------------------------------------------- 120 !! *** ROUTINE ldf_dyn_c2d *** 179 END SUBROUTINE ldf_dyn_c3d 180 181 182 SUBROUTINE ldf_dyn_c3d_orca( ld_print ) 183 !!---------------------------------------------------------------------- 184 !! *** ROUTINE ldf_dyn_c3d *** 185 !! 186 !! ** Purpose : ORCA R2 an R4 only 121 187 !! 122 !! **** W A R N I N G **** 123 !! 124 !! ORCA R2 and R4 configurations 125 !! 126 !! **** W A R N I N G **** 127 !! 128 !! ** Purpose : initializations of the lateral viscosity for orca R2 129 !! 130 !! ** Method : blah blah blah... 131 !! 188 !! ** Method : blah blah blah .... 132 189 !!---------------------------------------------------------------------- 133 190 !! * Modules used … … 135 192 136 193 !! * Arguments 137 LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout 138 139 !! * Local variables 140 INTEGER :: ji, jj, jn ! dummy loop indices 141 INTEGER :: inum ! temporary logical unit 194 LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout 195 196 !! * local variables 197 INTEGER :: ji, jj, jk, jn ! dummy loop indices 198 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 199 INTEGER :: inum ! temporary logical unit 142 200 INTEGER :: iim, ijm 143 201 INTEGER :: ifreq, il1, il2, ij, ii 144 INTEGER, DIMENSION(jpidta,jpidta) :: idata 145 INTEGER, DIMENSION(jpi ,jpj ) :: icof 146 147 REAL(wp) :: zahmeq, zcoft, zcoff, zmsk 202 INTEGER, DIMENSION(jpidta, jpjdta) :: idata 203 INTEGER, DIMENSION(jpi , jpj ) :: icof 204 205 REAL(wp) :: & 206 zahmeq, zcoff, zcoft, zmsk, & ! ??? 207 zemax, zemin, zeref, zahmm 208 REAL(wp), DIMENSION(jpi,jpj) :: zahm0 209 REAL(wp), DIMENSION(jpk) :: zcoef 148 210 149 211 CHARACTER (len=15) :: clexp … … 151 213 152 214 IF(lwp) WRITE(numout,*) 153 IF(lwp) WRITE(numout,*) ' inildf: 2deddy viscosity coefficient'154 IF(lwp) WRITE(numout,*) '~~~~~~ --'215 IF(lwp) WRITE(numout,*) 'ldfdyn_c3d_orca : 3D eddy viscosity coefficient' 216 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 155 217 IF(lwp) WRITE(numout,*) 156 IF(lwp) WRITE(numout,*) ' orca ocean model' 218 IF(lwp) WRITE(numout,*) ' orca R2 or R4 ocean model' 219 IF(lwp) WRITE(numout,*) ' reduced in the surface Eq. strip ' 157 220 IF(lwp) WRITE(numout,*) 158 221 159 #if defined key_antarctic160 # include "ldfdyn_antarctic.h90"161 #elif defined key_arctic162 # include "ldfdyn_arctic.h90"163 #else164 222 ! Read 2d integer array to specify western boundary increase in the 165 223 ! ===================== equatorial strip (20N-20S) defined at t-points 166 224 167 225 CALL ctlopn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', & 168 & 226 & 1, numout, lwp, 1 ) 169 227 REWIND inum 170 228 READ(inum,9101) clexp, iim, ijm … … 198 256 END DO 199 257 END DO 200 201 9101 FORMAT(1x,a15,2i8) 202 9201 FORMAT(3x,13(i3,12x)) 203 9202 FORMAT(i3,41i3) 204 205 206 ! Set ahm1 and ahm2 ( T- and F- points) (used for laplacian operator) 258 259 9101 FORMAT(1x,a15,2i8) 260 9201 FORMAT(3x,13(i3,12x)) 261 9202 FORMAT(i3,41i3) 262 263 ! Set ahm1 and ahm2 207 264 ! ================= 265 208 266 ! define ahm1 and ahm2 at the right grid point position 209 267 ! (USER: modify ahm1 and ahm2 following your desiderata) 210 268 ! biharmonic : ahm1 (ahm2) defined at u- (v-) point 269 ! harmonic : ahm1 (ahm2) defined at t- (f-) point 270 271 ! first level : as for 2D coefficients 211 272 212 273 ! Decrease ahm to zahmeq m2/s in the tropics … … 215 276 ! from 2.5 to 0 degre: ahm = constant 216 277 ! symmetric in the south hemisphere) 217 218 zahmeq = aht0 219 278 279 IF( jp_cfg == 4 ) THEN 280 zahmeq = 5.0 * aht0 281 zahmm = min( 160000.0, ahm0) 282 zemax = MAXVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 ) 283 zemin = MINVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 ) 284 zeref = MAXVAL ( e1t(:,:) * e2t(:,:), & 285 & tmask(:,:,1) .GE. 0.5 .AND. ABS(gphit(:,:)) .GT. 50. ) 286 287 DO jj = 1, jpj 288 DO ji = 1, jpi 289 zmsk = e1t(ji,jj) * e2t(ji,jj) 290 IF( abs(gphit(ji,jj)) .LE. 15 ) THEN 291 zahm0(ji,jj) = ahm0 292 ELSE 293 IF( zmsk .GE. zeref ) THEN 294 zahm0(ji,jj) = ahm0 295 ELSE 296 zahm0(ji,jj) = zahmm + (ahm0-zahmm)*(1.0 - & 297 & cos((rpi*0.5*(zmsk-zemin)/(zeref-zemin)))) 298 ENDIF 299 ENDIF 300 END DO 301 END DO 302 ENDIF 303 304 IF( jp_cfg == 2 ) THEN 305 zahmeq = aht0 306 zahmm = ahm0 307 zahm0(:,:) = ahm0 308 ENDIF 309 220 310 DO jj = 1, jpj 221 311 DO ji = 1, jpi 222 IF( ABS( gphif(ji,jj) ) >= 20.) THEN223 ahm2(ji,jj ) = ahm0224 ELSEIF( ABS( gphif(ji,jj) ) <= 2.5) THEN225 ahm2(ji,jj ) = zahmeq312 IF( ABS(gphif(ji,jj)) >= 20.) THEN 313 ahm2(ji,jj,1) = zahm0(ji,jj) 314 ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN 315 ahm2(ji,jj,1) = zahmeq 226 316 ELSE 227 ahm2(ji,jj ) = zahmeq + (ahm0-zahmeq)/2. &228 * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. /17.5 ) )317 ahm2(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. & 318 & *(1.-COS( rad*(ABS(gphif(ji,jj))-2.5)*180./17.5 ) ) 229 319 ENDIF 230 IF( ABS( gphit(ji,jj) ) >= 20.) THEN231 ahm1(ji,jj ) = ahm0232 ELSEIF( ABS( gphit(ji,jj) ) <= 2.5) THEN233 ahm1(ji,jj ) = zahmeq320 IF( ABS(gphit(ji,jj)) >= 20.) THEN 321 ahm1(ji,jj,1) = zahm0(ji,jj) 322 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 323 ahm1(ji,jj,1) = zahmeq 234 324 ELSE 235 ahm1(ji,jj ) = zahmeq + (ahm0-zahmeq)/2. &236 * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. /17.5 ) )325 ahm1(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. & 326 & *(1.-COS( rad*(ABS(gphit(ji,jj))-2.5)*180./17.5 ) ) 237 327 ENDIF 238 328 END DO 239 329 END DO 240 330 241 331 ! increase along western boundaries of equatorial strip 242 332 ! t-point 243 333 DO jj = 1, jpjm1 244 334 DO ji = 1, jpim1 245 zcoft = FLOAT( icof(ji,jj) ) / 100.246 ahm1(ji,jj ) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)335 zcoft = float( icof(ji,jj) ) / 100. 336 ahm1(ji,jj,1) = zcoft * zahm0(ji,jj) + (1.-zcoft) * ahm1(ji,jj,1) 247 337 END DO 248 338 END DO … … 258 348 / (zmsk * 100.) 259 349 ENDIF 260 ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj) 261 END DO 350 ahm2(ji,jj,1) = zcoff * zahm0(ji,jj) + (1.-zcoff) * ahm2(ji,jj,1) 351 END DO 352 END DO 353 354 ! other level: re-increase the coef in the deep ocean 355 356 #if defined key_orca_lev10 357 DO jk = 1, 210 358 zcoef(jk) = 1. 359 END DO 360 DO jk= 211, 230 361 zcoef(jk) = 1. + 0.1 * FLOAT(jk-210) 362 END DO 363 DO jk= 231, 260 364 zcoef(jk) = 3. + 0.2 * FLOAT(jk-230) 365 END DO 366 DO jk= 261, 270 367 zcoef(jk) = 9. + 0.1 * FLOAT(jk-260) 368 END DO 369 DO jk= 271, jpk 370 zcoef(jk) = 10. 371 END DO 372 DO jk= 1, jpk 373 IF(lwp) WRITE(numout,*) 'k= ',jk, 'cof ', zcoef(jk) 374 END DO 375 #else 376 DO jk = 1, 21 377 zcoef(jk) = 1. 378 END DO 379 zcoef(22) = 2. 380 zcoef(23) = 3. 381 zcoef(24) = 5. 382 zcoef(25) = 7. 383 zcoef(26) = 9. 384 DO jk = 27, jpk 385 zcoef(jk) = 10. 262 386 END DO 263 387 #endif 388 389 DO jk = 2, jpk 390 ahm1(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm1(:,:,1) ) 391 ahm2(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm2(:,:,1) ) 392 END DO 393 394 IF( jp_cfg == 4 ) THEN 395 ! Limit AHM in Gibraltar strait 396 ij0 = 50 ; ij1 = 53 397 ii0 = 69 ; ii1 = 71 398 DO jk = 1, jpk 399 ahm1(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) = min( zahmm, ahm1 (mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) ) 400 ahm2(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) = min( zahmm, ahm2 (mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) ) 401 END DO 402 ENDIF 264 403 265 404 ! Lateral boundary conditions on ( ahm1, ahm2 ) … … 269 408 270 409 ! Control print 410 411 IF(lwp) THEN 412 WRITE(numout,*) 413 WRITE(numout,*) ' 3D ahm1 array (k=1)' 414 CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout ) 415 WRITE(numout,*) 416 WRITE(numout,*) ' 3D ahm2 array (k=1)' 417 CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout ) 418 WRITE(numout,*) 419 WRITE(numout,*) ' 3D ahm2 array (k=jpk)' 420 CALL prihre( ahm2(:,:,jpk), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout ) 421 ENDIF 422 423 424 ! Set ahm3 and ahm4 425 ! ================= 426 427 ! define ahm3 and ahm4 at the right grid point position 428 ! initialization to a constant value 429 ! (USER: modify ahm3 and ahm4 following your desiderata) 430 ! harmonic isopycnal or geopotential: 431 ! ahm3 (ahm4) defined at u- (v-) point 432 DO jk = 1, jpk 433 DO jj = 2, jpj 434 DO ji = 2, jpi 435 ahm3(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji ,jj-1,jk) ) 436 ahm4(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji-1,jj ,jk) ) 437 END DO 438 END DO 439 END DO 440 ahm3 ( :, 1, :) = ahm3 ( :, 2, :) 441 ahm4 ( :, 1, :) = ahm4 ( :, 2, :) 442 443 ! Lateral boundary conditions on ( ahm3, ahm4 ) 444 ! ============== 445 CALL lbc_lnk( ahm3, 'U', 1. ) ! U-point, unchanged sign 446 CALL lbc_lnk( ahm4, 'V', 1. ) ! V-point, unchanged sign 447 448 ! Control print 449 271 450 IF( lwp .AND. ld_print ) THEN 272 451 WRITE(numout,*) 273 WRITE(numout,*) ' inildf: 2D ahm1 array'274 CALL prihre(ahm 1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)275 WRITE(numout,*) 276 WRITE(numout,*) ' inildf: 2D ahm2 array'277 CALL prihre(ahm 2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)278 ENDIF 279 280 END SUBROUTINE ldf_dyn_c 2d_orca452 WRITE(numout,*) ' ahm3 array level 1' 453 CALL prihre(ahm3(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 454 WRITE(numout,*) 455 WRITE(numout,*) ' ahm4 array level 1' 456 CALL prihre(ahm4(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 457 ENDIF 458 459 END SUBROUTINE ldf_dyn_c3d_orca
Note: See TracChangeset
for help on using the changeset viewer.
