- Timestamp:
- 2017-12-13T15:37:41+01:00 (6 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
r8995 r9017 1 MODULE usrdef_closea1 MODULE closea 2 2 !!====================================================================== 3 3 !! *** MODULE usrdef_closea *** 4 !!5 !! === ORCA configuration ===6 !! (2, 1 and 1/4 degrees)7 4 !! 8 5 !! User define : specific treatments associated with closed seas … … 13 10 !! 3.4 ! 2014-12 (P.G. Fogli) sbc_clo bug fix & mpp reproducibility 14 11 !! 4.0 ! 2016-06 (G. Madec) move to usrdef_closea, remove clo_ups 12 !! 4.0 ! 2017-12 (D. Storkey) new formulation based on masks read from file 15 13 !!---------------------------------------------------------------------- 16 14 … … 25 23 USE phycst ! physical constants 26 24 USE sbc_oce ! ocean surface boundary conditions 25 USE iom ! I/O routines 27 26 ! 28 27 USE in_out_manager ! I/O manager 29 USE lib_fortran, ONLY: glob_sum , DDPDD28 USE lib_fortran, ONLY: glob_sum 30 29 USE lbclnk ! lateral boundary condition - MPP exchanges 31 30 USE lib_mpp ! MPP library 32 31 USE timing 32 USE wrk_nemo ! Memory allocation 33 33 34 34 IMPLICIT NONE … … 36 36 37 37 PUBLIC dom_clo ! called by domain module 38 PUBLIC sbc_clo ! called by s tepmodule38 PUBLIC sbc_clo ! called by sbcmod module 39 39 PUBLIC clo_rnf ! called by sbcrnf module 40 PUBLIC clo_bat ! called in domzgr module 41 42 INTEGER, PUBLIC, PARAMETER :: jpncs = 4 !: number of closed sea 43 INTEGER, PUBLIC, DIMENSION(jpncs) :: ncstt !: Type of closed sea 44 INTEGER, PUBLIC, DIMENSION(jpncs) :: ncsi1, ncsj1 !: south-west closed sea limits (i,j) 45 INTEGER, PUBLIC, DIMENSION(jpncs) :: ncsi2, ncsj2 !: north-east closed sea limits (i,j) 46 INTEGER, PUBLIC, DIMENSION(jpncs) :: ncsnr !: number of point where run-off pours 47 INTEGER, PUBLIC, DIMENSION(jpncs,4) :: ncsir, ncsjr !: Location of runoff 48 49 REAL(wp), DIMENSION (jpncs+1) :: surf ! closed sea surface 40 PUBLIC clo_bat ! called in domain module 41 42 INTEGER, PUBLIC :: jncs !: number of closed seas (inferred from closea_mask field) 43 INTEGER, PUBLIC :: jncsr !: number of closed seas rnf mappings (inferred from closea_mask_rnf field) 44 INTEGER, PUBLIC :: jncse !: number of closed seas empmr mappings (inferred from closea_mask_empmr field) 45 46 INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: closea_mask !: mask of integers defining closed seas 47 INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: closea_mask_rnf !: mask of integers defining closed seas rnf mappings 48 INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: closea_mask_empmr !: mask of integers defining closed seas empmr mappings 49 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: surf !: closed sea surface areas 50 !: (and residual global surface area) 51 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: surfr !: closed sea target rnf surface areas 52 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: surfe !: closed sea target empmr surface areas 50 53 51 54 !! * Substitutions … … 58 61 CONTAINS 59 62 60 SUBROUTINE dom_clo( cd_cfg, kcfg)63 SUBROUTINE dom_clo( k_bot ) 61 64 !!--------------------------------------------------------------------- 62 65 !! *** ROUTINE dom_clo *** … … 67 70 !! just the thermodynamic processes are applied. 68 71 !! 69 !! ** Action : ncsi1(), ncsj1() : south-west closed sea limits (i,j) 70 !! ncsi2(), ncsj2() : north-east Closed sea limits (i,j) 71 !! ncsir(), ncsjr() : Location of runoff 72 !! ncsnr : number of point where run-off pours 73 !! ncstt : Type of closed sea 74 !! =0 spread over the world ocean 75 !! =2 put at location runoff 76 !!---------------------------------------------------------------------- 77 CHARACTER(len=*), INTENT(in ) :: cd_cfg ! configuration name 78 INTEGER , INTENT(in ) :: kcfg ! configuration identifier 79 ! 80 INTEGER :: jc ! dummy loop indices 81 INTEGER :: isrow ! local index 72 !! ** Action : Read closea mask fields from closea_masks.nc and infer 73 !! number of closed seas from closea_mask field. 74 !! closea_mask : integer values defining closed seas (or groups of closed seas) 75 !! closea_mask_rnf : integer values defining mappings from closed seas or groups of 76 !! closed seas to a runoff area for downwards flux only. 77 !! closea_mask_empmr : integer values defining mappings from closed seas or groups of 78 !! closed seas to a runoff area for net fluxes. 79 !! 80 !! Python code to generate the closea_masks.nc file from the old-style indices 81 !! definitions is available at TOOLS/DOMAINcfg/make_closea_masks.py 82 !!---------------------------------------------------------------------- 83 INTEGER, DIMENSION(:,:), INTENT(in) :: k_bot ! ocean last level index (for masking input fields) 84 INTEGER :: inum ! input file identifier 85 INTEGER :: ierr ! error code 86 INTEGER :: id ! netcdf variable ID 87 REAL(wp), POINTER, DIMENSION(:,:) :: zdata_in ! temporary real array for input 82 88 !!---------------------------------------------------------------------- 83 89 ! … … 86 92 IF(lwp) WRITE(numout,*)'~~~~~~~' 87 93 ! 88 ! initial values 89 ncsnr(:) = 1 ; ncsi1(:) = 1 ; ncsi2(:) = 1 ; ncsir(:,:) = 1 90 ncstt(:) = 0 ; ncsj1(:) = 1 ; ncsj2(:) = 1 ; ncsjr(:,:) = 1 91 ! 92 ! set the closed seas (in data domain indices) 93 ! ------------------- 94 ! 95 IF( cd_cfg == "orca" ) THEN !== ORCA configuration ==! 96 ! 97 SELECT CASE ( kcfg ) 98 ! ! ======================= 99 CASE ( 1 ) ! ORCA_R1 configuration 100 ! ! ======================= 101 IF(lwp) WRITE(numout,*)' ORCA_R1 closed seas : only the Caspian Sea' 102 ! This dirty section will be suppressed by simplification process: 103 ! all this will come back in input files 104 ! Currently these hard-wired indices relate to configuration with 105 ! extend grid (jpjglo=332) 106 isrow = 332 - jpjglo 107 ! 108 ncsnr(1) = 1 ; ncstt(1) = 0 ! Caspian Sea (spread over the globe) 109 ncsi1(1) = 332 ; ncsj1(1) = 243 - isrow 110 ncsi2(1) = 344 ; ncsj2(1) = 275 - isrow 111 ncsir(1,1) = 1 ; ncsjr(1,1) = 1 112 ! 113 ! ! ======================= 114 CASE ( 2 ) ! ORCA_R2 configuration 115 ! ! ======================= 116 IF(lwp) WRITE(numout,*)' ORCA_R2 closed seas and lakes : ' 117 ! ! Caspian Sea 118 IF(lwp) WRITE(numout,*)' Caspian Sea ' 119 ncsnr(1) = 1 ; ncstt(1) = 0 ! spread over the globe 120 ncsi1(1) = 11 ; ncsj1(1) = 103 121 ncsi2(1) = 17 ; ncsj2(1) = 112 122 ncsir(1,1) = 1 ; ncsjr(1,1) = 1 123 ! ! Great North American Lakes 124 IF(lwp) WRITE(numout,*)' Great North American Lakes ' 125 ncsnr(2) = 1 ; ncstt(2) = 2 ! put at St Laurent mouth 126 ncsi1(2) = 97 ; ncsj1(2) = 107 127 ncsi2(2) = 103 ; ncsj2(2) = 111 128 ncsir(2,1) = 110 ; ncsjr(2,1) = 111 129 ! ! Black Sea (crossed by the cyclic boundary condition) 130 IF(lwp) WRITE(numout,*)' Black Sea ' 131 ncsnr(3:4) = 4 ; ncstt(3:4) = 2 ! put in Med Sea (north of Aegean Sea) 132 ncsir(3:4,1) = 171; ncsjr(3:4,1) = 106 ! 133 ncsir(3:4,2) = 170; ncsjr(3:4,2) = 106 134 ncsir(3:4,3) = 171; ncsjr(3:4,3) = 105 135 ncsir(3:4,4) = 170; ncsjr(3:4,4) = 105 136 ncsi1(3) = 174 ; ncsj1(3) = 107 ! 1 : west part of the Black Sea 137 ncsi2(3) = 181 ; ncsj2(3) = 112 ! (ie west of the cyclic b.c.) 138 ncsi1(4) = 2 ; ncsj1(4) = 107 ! 2 : east part of the Black Sea 139 ncsi2(4) = 6 ; ncsj2(4) = 112 ! (ie east of the cyclic b.c.) 140 ! 141 ! ! ========================= 142 CASE ( 025 ) ! ORCA_R025 configuration 143 ! ! ========================= 144 IF(lwp) WRITE(numout,*)' ORCA_R025 closed seas : ' 145 ! ! Caspian Sea 146 IF(lwp) WRITE(numout,*)' Caspian Sea ' 147 ncsnr(1) = 1 ; ncstt(1) = 0 ! Caspian + Aral sea 148 ncsi1(1) = 1330 ; ncsj1(1) = 645 149 ncsi2(1) = 1400 ; ncsj2(1) = 795 150 ncsir(1,1) = 1 ; ncsjr(1,1) = 1 151 ! 152 IF(lwp) WRITE(numout,*)' Azov Sea ' 153 ncsnr(2) = 1 ; ncstt(2) = 0 ! Azov Sea 154 ncsi1(2) = 1284 ; ncsj1(2) = 722 155 ncsi2(2) = 1304 ; ncsj2(2) = 747 156 ncsir(2,1) = 1 ; ncsjr(2,1) = 1 157 ! 158 END SELECT 159 ! 160 ELSE !== No closed sea in the configuration ==! 161 ! 162 IF(lwp) WRITE(numout,*)' No closed seas or lakes in the configuration ' 163 ! 94 CALL wrk_alloc( jpi,jpj, zdata_in) 95 ! 96 ! read the closed seas masks 97 ! -------------------------- 98 ! 99 ALLOCATE( closea_mask(jpi,jpj) , STAT=ierr ) 100 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask array') 101 ALLOCATE( closea_mask_rnf(jpi,jpj) , STAT=ierr ) 102 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_rnf array') 103 ALLOCATE( closea_mask_empmr(jpi,jpj) , STAT=ierr ) 104 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_empmr array') 105 106 closea_mask(:,:) = 0 107 closea_mask_rnf(:,:) = 0 108 closea_mask_empmr(:,:) = 0 109 110 CALL iom_open( 'closea_masks.nc', inum ) 111 112 zdata_in(:,:) = 0.0 113 CALL iom_get ( inum, jpdom_data, 'closea_mask', zdata_in ) 114 closea_mask(:,:) = NINT(zdata_in(:,:)) 115 ! necessary because fill values in input fields can confuse things 116 ! we can't multiply by tmask here because it isn't defined yet 117 WHERE( k_bot(:,:) == 0 ) closea_mask(:,:) = 0 118 119 id = iom_varid(inum, 'closea_mask_rnf', ldstop = .false.) 120 IF(lwp) WRITE(numout,*) 'closea_mask_rnf, id : ',id 121 IF( id > 0 ) THEN 122 CALL iom_get ( inum, jpdom_data, 'closea_mask_rnf', zdata_in ) 123 closea_mask_rnf(:,:) = NINT(zdata_in(:,:)) 124 WHERE( k_bot(:,:) == 0 ) closea_mask_rnf(:,:) = 0 164 125 ENDIF 165 166 ! convert the position in local domain indices 167 ! -------------------------------------------- 168 DO jc = 1, jpncs 169 ncsi1(jc) = mi0( ncsi1(jc) ) 170 ncsj1(jc) = mj0( ncsj1(jc) ) 171 ! 172 ncsi2(jc) = mi1( ncsi2(jc) ) 173 ncsj2(jc) = mj1( ncsj2(jc) ) 174 END DO 126 127 id = iom_varid(inum, 'closea_mask_empmr', ldstop = .false.) 128 IF( id > 0 ) THEN 129 CALL iom_get ( inum, jpdom_data, 'closea_mask_empmr', zdata_in ) 130 closea_mask_empmr(:,:) = NINT(zdata_in(:,:)) 131 WHERE( k_bot(:,:) == 0 ) closea_mask_empmr(:,:) = 0 132 ENDIF 133 134 CALL iom_close( inum ) 135 136 ! number of closed seas = global maximum value in closea_mask field 137 jncs = maxval(closea_mask(:,:)) 138 IF( lk_mpp ) CALL mpp_max(jncs) 139 IF( lwp ) WRITE(numout,*) 'Number of closed seas : ',jncs 140 ! 141 ! number of closed seas rnf mappings = global maximum in closea_mask_rnf field 142 jncsr = maxval(closea_mask_rnf(:,:)) 143 IF( lk_mpp ) CALL mpp_max(jncsr) 144 IF( lwp ) WRITE(numout,*) 'Number of closed seas rnf mappings : ',jncsr 145 ! 146 ! number of closed seas empmr mappings = global maximum value in closea_mask_empmr field 147 jncse = maxval(closea_mask_empmr(:,:)) 148 IF( lk_mpp ) CALL mpp_max(jncse) 149 IF( lwp ) WRITE(numout,*) 'Number of closed seas empmr mappings : ',jncse 150 ! 151 CALL wrk_dealloc(jpi,jpj, zdata_in) 175 152 ! 176 153 END SUBROUTINE dom_clo 177 154 178 155 179 SUBROUTINE sbc_clo( kt , cd_cfg, kcfg)156 SUBROUTINE sbc_clo( kt ) 180 157 !!--------------------------------------------------------------------- 181 158 !! *** ROUTINE sbc_clo *** … … 190 167 !!---------------------------------------------------------------------- 191 168 INTEGER , INTENT(in ) :: kt ! ocean model time step 192 CHARACTER(len=*), INTENT(in ) :: cd_cfg ! configuration name 193 INTEGER , INTENT(in ) :: kcfg ! configuration identifier 194 ! 195 INTEGER :: ji, jj, jc, jn ! dummy loop indices 169 ! 170 INTEGER :: ierr 171 INTEGER :: jc, jcr, jce ! dummy loop indices 196 172 REAL(wp), PARAMETER :: rsmall = 1.e-20_wp ! Closed sea correction epsilon 197 REAL(wp) :: zze2, ztmp, zcorr ! 198 REAL(wp) :: zcoef, zcoef1 ! 199 COMPLEX(wp) :: ctmp 200 REAL(wp), DIMENSION(jpncs) :: zfwf ! 1D workspace 173 REAL(wp) :: zfwf_total, zcoef, zcoef1 ! 174 REAL(wp), POINTER, DIMENSION(:) :: zfwf, zfwfr, zfwfe ! freshwater fluxes over closed seas 175 REAL(wp), POINTER, DIMENSION(:,:) :: ztmp2d ! 2D workspace 201 176 !!---------------------------------------------------------------------- 202 177 ! 203 178 IF( nn_timing == 1 ) CALL timing_start('sbc_clo') 204 ! !------------------! 179 ! 180 CALL wrk_alloc( jncs, zfwf ) 181 ! jncsr and jncse can be zero so add 1 to avoid allocating zero-length arrays 182 CALL wrk_alloc( jncsr+1, zfwfr ) 183 CALL wrk_alloc( jncse+1, zfwfe ) 184 CALL wrk_alloc( jpi,jpj, ztmp2d ) 185 ! !------------------! 205 186 IF( kt == nit000 ) THEN ! Initialisation ! 206 187 ! !------------------! … … 209 190 IF(lwp) WRITE(numout,*)'~~~~~~~' 210 191 211 surf(:) = 0._wp 212 ! 213 surf(jpncs+1) = glob_sum( e1e2t(:,:) ) ! surface of the global ocean 214 ! 215 ! ! surface of closed seas 216 DO jc = 1, jpncs 217 ctmp = CMPLX( 0.e0, 0.e0, wp ) 218 DO jj = ncsj1(jc), ncsj2(jc) 219 DO ji = ncsi1(jc), ncsi2(jc) 220 ztmp = e1e2t(ji,jj) * tmask_i(ji,jj) 221 CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 222 END DO 223 END DO 224 IF( lk_mpp ) CALL mpp_sum( ctmp ) 225 surf(jc) = REAL(ctmp,wp) 192 ALLOCATE( surf(jncs+1) , STAT=ierr ) 193 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surf array') 194 surf(:) = 0.e0_wp 195 ! 196 ! jncsr can be zero so add 1 to avoid allocating zero-length array 197 ALLOCATE( surfr(jncsr+1) , STAT=ierr ) 198 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surfr array') 199 surfr(:) = 0.e0_wp 200 ! 201 ! jncse can be zero so add 1 to avoid allocating zero-length array 202 ALLOCATE( surfe(jncse+1) , STAT=ierr ) 203 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surfe array') 204 surfe(:) = 0.e0_wp 205 ! 206 surf(jncs+1) = glob_sum( e1e2t(:,:) ) ! surface of the global ocean 207 ! 208 ! ! surface areas of closed seas 209 DO jc = 1, jncs 210 ztmp2d(:,:) = 0.e0_wp 211 WHERE( closea_mask(:,:) == jc ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:) 212 surf(jc) = glob_sum( ztmp2d(:,:) ) 226 213 END DO 227 228 IF(lwp) WRITE(numout,*)' Closed sea surfaces' 229 DO jc = 1, jpncs 230 IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc) 214 ! 215 ! jncs+1 : surface area of global ocean, closed seas excluded 216 surf(jncs+1) = surf(jncs+1) - SUM(surf(1:jncs)) 217 ! 218 ! ! surface areas of rnf target areas 219 DO jcr = 1, jncsr 220 ztmp2d(:,:) = 0.e0_wp 221 WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:) 222 surfr(jcr) = glob_sum( ztmp2d(:,:) ) 231 223 END DO 232 233 ! jpncs+1 : surface of sea, closed seas excluded 234 DO jc = 1, jpncs 235 surf(jpncs+1) = surf(jpncs+1) - surf(jc) 236 END DO 237 ! 224 ! 225 ! ! surface areas of empmr target areas 226 DO jce = 1, jncse 227 ztmp2d(:,:) = 0.e0_wp 228 WHERE( closea_mask_empmr(:,:) == jcr .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:) 229 surfe(jcr) = glob_sum( ztmp2d(:,:) ) 230 END DO 231 ! 232 IF(lwp) WRITE(numout,*)' Closed sea surface areas (km2)' 233 DO jc = 1, jncs 234 IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jc, surf(jc) * 1.0e-6 235 END DO 236 IF(lwp) WRITE(numout,FMT='(A,ES12.2)') 'Global surface area excluding closed seas (km2): ', surf(jncs+1) * 1.0e-6 237 ! 238 IF(jncsr > 0) THEN 239 IF(lwp) WRITE(numout,*)' Closed sea target rnf surface areas (km2)' 240 DO jcr = 1, jncsr 241 IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jcr, surfr(jcr) * 1.0e-6 242 END DO 243 ENDIF 244 ! 245 IF(jncse > 0) THEN 246 IF(lwp) WRITE(numout,*)' Closed sea target empmr surface areas (km2)' 247 DO jce = 1, jncse 248 IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jce, surfe(jce) * 1.0e-6 249 END DO 250 ENDIF 238 251 ENDIF 239 ! !--------------------! 240 ! ! update emp ! 241 zfwf = 0.e0_wp !--------------------! 242 DO jc = 1, jpncs 243 ctmp = CMPLX( 0.e0, 0.e0, wp ) 244 DO jj = ncsj1(jc), ncsj2(jc) 245 DO ji = ncsi1(jc), ncsi2(jc) 246 ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 247 CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 248 END DO 249 END DO 250 IF( lk_mpp ) CALL mpp_sum( ctmp ) 251 zfwf(jc) = REAL(ctmp,wp) 252 ! 253 ! !--------------------! 254 ! ! update emp ! 255 ! !--------------------! 256 257 zfwf_total = 0._wp 258 259 ! 260 ! 1. Work out total freshwater fluxes over closed seas from EMP - RNF. 261 ! 262 zfwf(:) = 0.e0_wp 263 DO jc = 1, jncs 264 ztmp2d(:,:) = 0.e0_wp 265 WHERE( closea_mask(:,:) == jc ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) 266 zfwf(jc) = glob_sum( ztmp2d(:,:) ) 252 267 END DO 253 254 IF( cd_cfg == "orca" .AND. kcfg == 2 ) THEN ! Black Sea case for ORCA_R2 configuration 255 zze2 = ( zfwf(3) + zfwf(4) ) * 0.5_wp 256 zfwf(3) = zze2 257 zfwf(4) = zze2 268 zfwf_total = SUM(zfwf) 269 270 ! 271 ! 2. Work out total FW fluxes over rnf source areas and add to rnf target areas. If jncsr is zero does not loop. 272 ! Where zfwf is negative add flux at specified runoff points and subtract from fluxes for global redistribution. 273 ! Where positive leave in global redistribution total. 274 ! 275 zfwfr(:) = 0.e0_wp 276 DO jcr = 1, jncsr 277 ! 278 ztmp2d(:,:) = 0.e0_wp 279 WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) > 0 ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) 280 zfwfr(jcr) = glob_sum( ztmp2d(:,:) ) 281 IF(lwp) WRITE(numout,*) 'closea runoff output: jcr, zfwfr(jcr), ABS(zfwfr(jcr)/surf(jncs+1) : ',jcr,zfwfr(jcr),ABS(zfwfr(jcr)/surf(jncs+1)) 282 ! 283 ! The following if avoids the redistribution of the round off 284 IF ( ABS(zfwfr(jcr) / surf(jncs+1) ) > rsmall) THEN 285 ! 286 ! Add residuals to target runoff points if negative and subtract from total to be added globally 287 IF(lwp) WRITE(numout,*) 'closea runoff output: passed roundoff test!' 288 IF( zfwfr(jcr) < 0.0 ) THEN 289 zfwf_total = zfwf_total - zfwfr(jcr) 290 zcoef = zfwfr(jcr) / surfr(jcr) 291 zcoef1 = rcp * zcoef 292 IF(lwp) WRITE(numout,*) 'closea runoff output: jcr, zcoef: ',jcr, zcoef 293 WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0.0) 294 emp(:,:) = emp(:,:) - zcoef 295 qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:) 296 ENDWHERE 297 ENDIF 298 ! 299 ENDIF 300 END DO 301 302 ! 303 ! 3. Work out total fluxes over empmr source areas and add to empmr target areas. If jncse is zero does not loop. 304 ! 305 zfwfe(:) = 0.e0_wp 306 DO jce = 1, jncse 307 ! 308 ztmp2d(:,:) = 0.e0_wp 309 WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) > 0 ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) 310 zfwfe(jce) = glob_sum( ztmp2d(:,:) ) 311 ! 312 ! The following if avoids the redistribution of the round off 313 IF ( ABS(zfwfe(jce) / surf(jncs+1) ) > rsmall) THEN 314 ! 315 ! Add residuals to runoff points and subtract from total to be added globally 316 zfwf_total = zfwf_total - zfwfe(jce) 317 zcoef = zfwfe(jce) / surfe(jce) 318 zcoef1 = rcp * zcoef 319 WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) == 0.0) 320 emp(:,:) = emp(:,:) - zcoef 321 qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:) 322 ENDWHERE 323 ! 324 ENDIF 325 END DO 326 327 ! 328 ! 4. Spread residual flux over global ocean. 329 ! 330 ! The following if avoids the redistribution of the round off 331 IF ( ABS(zfwf_total / surf(jncs+1) ) > rsmall) THEN 332 zcoef = zfwf_total / surf(jncs+1) 333 zcoef1 = rcp * zcoef 334 IF(lwp) WRITE(numout,*) 'closea global addition: zfwf_total, zcoef: ', zfwf_total, zcoef 335 WHERE( closea_mask(:,:) == 0 ) 336 emp(:,:) = emp(:,:) + zcoef 337 qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 338 ENDWHERE 258 339 ENDIF 259 340 260 zcorr = 0._wp261 262 DO jc = 1, jpncs263 !341 ! 342 ! 5. Subtract area means from emp (and qns) over closed seas to give zero mean FW flux over each sea. 343 ! 344 DO jc = 1, jncs 264 345 ! The following if avoids the redistribution of the round off 265 IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN 266 ! 267 IF( ncstt(jc) == 0 ) THEN ! water/evap excess is shared by all open ocean 268 zcoef = zfwf(jc) / surf(jpncs+1) 269 zcoef1 = rcp * zcoef 270 emp(:,:) = emp(:,:) + zcoef 271 qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 272 ! accumulate closed seas correction 273 zcorr = zcorr + zcoef 274 ! 275 ELSEIF( ncstt(jc) == 1 ) THEN ! Excess water in open sea, at outflow location, excess evap shared 276 IF ( zfwf(jc) <= 0.e0_wp ) THEN 277 DO jn = 1, ncsnr(jc) 278 ji = mi0(ncsir(jc,jn)) 279 jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 280 IF ( ji > 1 .AND. ji < jpi & 281 .AND. jj > 1 .AND. jj < jpj ) THEN 282 zcoef = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 283 zcoef1 = rcp * zcoef 284 emp(ji,jj) = emp(ji,jj) + zcoef 285 qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj) 286 ENDIF 287 END DO 288 ELSE 289 zcoef = zfwf(jc) / surf(jpncs+1) 290 zcoef1 = rcp * zcoef 291 emp(:,:) = emp(:,:) + zcoef 292 qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 293 ! accumulate closed seas correction 294 zcorr = zcorr + zcoef 295 ENDIF 296 ELSEIF( ncstt(jc) == 2 ) THEN ! Excess e-p-r (either sign) goes to open ocean, at outflow location 297 DO jn = 1, ncsnr(jc) 298 ji = mi0(ncsir(jc,jn)) 299 jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 300 IF( ji > 1 .AND. ji < jpi & 301 .AND. jj > 1 .AND. jj < jpj ) THEN 302 zcoef = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 303 zcoef1 = rcp * zcoef 304 emp(ji,jj) = emp(ji,jj) + zcoef 305 qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj) 306 ENDIF 307 END DO 308 ENDIF 309 ! 310 DO jj = ncsj1(jc), ncsj2(jc) 311 DO ji = ncsi1(jc), ncsi2(jc) 312 zcoef = zfwf(jc) / surf(jc) 313 zcoef1 = rcp * zcoef 314 emp(ji,jj) = emp(ji,jj) - zcoef 315 qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj) 316 END DO 317 END DO 318 ! 319 END IF 320 END DO 321 322 IF ( ABS(zcorr) > rsmall ) THEN ! remove the global correction from the closed seas 323 DO jc = 1, jpncs ! only if it is large enough 324 DO jj = ncsj1(jc), ncsj2(jc) 325 DO ji = ncsi1(jc), ncsi2(jc) 326 emp(ji,jj) = emp(ji,jj) - zcorr 327 qns(ji,jj) = qns(ji,jj) + rcp * zcorr * sst_m(ji,jj) 328 END DO 329 END DO 330 END DO 331 ENDIF 346 IF ( ABS(zfwf(jc) / surf(jncs+1) ) > rsmall) THEN 347 ! 348 ! Subtract residuals from fluxes over closed sea 349 zcoef = zfwf(jc) / surf(jc) 350 zcoef1 = rcp * zcoef 351 IF(lwp) WRITE(numout,*) 'closea subtract mean: jc, zfwf(jc), zcoef: ',jc, zfwf(jc), zcoef 352 WHERE( closea_mask(:,:) == jc ) 353 emp(:,:) = emp(:,:) - zcoef 354 qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:) 355 ENDWHERE 356 ! 357 ENDIF 358 END DO 332 359 ! 333 360 emp (:,:) = emp (:,:) * tmask(:,:,1) … … 335 362 CALL lbc_lnk( emp , 'T', 1._wp ) 336 363 ! 364 CALL wrk_dealloc( jncs , zfwf ) 365 CALL wrk_dealloc( jncsr+1 , zfwfr ) 366 CALL wrk_dealloc( jncse+1 , zfwfe ) 367 CALL wrk_dealloc( jpi,jpj, ztmp2d ) 368 ! 337 369 IF( nn_timing == 1 ) CALL timing_stop('sbc_clo') 338 370 ! 339 371 END SUBROUTINE sbc_clo 340 341 372 342 373 SUBROUTINE clo_rnf( p_rnfmsk ) … … 353 384 !!---------------------------------------------------------------------- 354 385 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: p_rnfmsk ! river runoff mask (rnfmsk array) 355 ! 356 INTEGER :: jc, jn, ji, jj ! dummy loop indices 357 !!---------------------------------------------------------------------- 358 ! 359 DO jc = 1, jpncs 360 IF( ncstt(jc) >= 1 ) THEN ! runoff mask set to 1 at closed sea outflows 361 DO jn = 1, 4 362 DO jj = mj0( ncsjr(jc,jn) ), mj1( ncsjr(jc,jn) ) 363 DO ji = mi0( ncsir(jc,jn) ), mi1( ncsir(jc,jn) ) 364 p_rnfmsk(ji,jj) = MAX( p_rnfmsk(ji,jj), 1.0_wp ) 365 END DO 366 END DO 367 END DO 368 ENDIF 369 END DO 386 !!---------------------------------------------------------------------- 387 ! 388 WHERE( ( closea_mask_rnf(:,:) > 0 .or. closea_mask_empmr(:,:) > 0 ) .and. closea_mask(:,:) == 0 ) 389 p_rnfmsk(:,:) = MAX( p_rnfmsk(:,:), 1.0_wp ) 390 ENDWHERE 370 391 ! 371 392 END SUBROUTINE clo_rnf … … 383 404 !!---------------------------------------------------------------------- 384 405 INTEGER, DIMENSION(:,:), INTENT(inout) :: k_top, k_bot ! ocean first and last level indices 385 ! 386 INTEGER :: jc, ji, jj ! dummy loop indices 387 !!---------------------------------------------------------------------- 388 ! 389 DO jc = 1, jpncs 390 DO jj = ncsj1(jc), ncsj2(jc) 391 DO ji = ncsi1(jc), ncsi2(jc) 392 k_top(ji,jj) = 0 393 k_bot(ji,jj) = 0 394 END DO 395 END DO 396 END DO 397 ! 406 !!---------------------------------------------------------------------- 407 ! 408 WHERE( closea_mask(:,:) > 0 ) 409 k_top(:,:) = 0 410 k_bot(:,:) = 0 411 ENDWHERE 412 ! 398 413 END SUBROUTINE clo_bat 399 414 400 415 !!====================================================================== 401 END MODULE usrdef_closea402 416 END MODULE closea 417
Note: See TracChangeset
for help on using the changeset viewer.
