- Timestamp:
- 2018-07-26T09:50:51+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90
r9939 r10001 14 14 USE dom_oce ! ocean space and time domain 15 15 USE domwri ! ocean space and time domain 16 USE domvvl , ONLY : dom_vvl_interpol17 16 USE phycst ! physical constants 18 17 USE sbc_oce ! surface boundary condition variables … … 137 136 REAL(wp), DIMENSION(jpi,jpj) :: zdmask , zsmask0, zucorr, zbub, zbun, zssh0, zhu1, zde3t 138 137 REAL(wp), DIMENSION(jpi,jpj) :: zdsmask, zsmask1, zvcorr, zbvb, zbvn, zssh1, zhv1 139 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask0, zwmaskn , ztrp138 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask0, zwmaskn 140 139 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask1, zwmaskb, ztmp3d 141 140 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0 141 REAL(wp), DIMENSION(jpi,jpj) :: z_ssh_h0, zsshu, zsshv, zsshf 142 142 !!---------------------------------------------------------------------- 143 143 ! … … 163 163 DO jj = 2,jpj-1 164 164 DO ji = fs_2, fs_jpim1 ! vector opt. 165 jip1=ji+1; jim1=ji-1; 166 jjp1=jj+1; jjm1=jj-1; 167 summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1)) 165 summsk = zsmask0(ji+1,jj)+zsmask0(ji-1,jj)+zsmask0(ji,jj+1)+zsmask0(ji,jj-1) 168 166 IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN 169 sshn(ji,jj)=( zssh0(ji p1,jj)*zsmask0(jip1,jj) &170 & + zssh0(ji m1,jj)*zsmask0(jim1,jj) &171 & + zssh0(ji,jj p1)*zsmask0(ji,jjp1) &172 & + zssh0(ji,jj m1)*zsmask0(ji,jjm1))/summsk173 zsmask1(ji,jj) =1._wp167 sshn(ji,jj)=( zssh0(ji+1,jj)*zsmask0(ji+1,jj) & 168 & + zssh0(ji-1,jj)*zsmask0(ji-1,jj) & 169 & + zssh0(ji,jj+1)*zsmask0(ji,jj+1) & 170 & + zssh0(ji,jj-1)*zsmask0(ji,jj+1)) / summsk 171 zsmask1(ji,jj) = 1._wp 174 172 ENDIF 175 173 END DO … … 185 183 IF ( .NOT.ln_linssh ) THEN 186 184 ! Reconstruction of all vertical scale factors at now time steps 187 ! ============================================================================= 185 ! ====================================================================== 186 187 !!gm Question : bug ???? 188 ! at this stage is ht_0, r1_ht0 already computed ???? 189 ! I have the feeling that this should be done in a different manner... 190 ! that is iscpl_rst_interpol should be call directly in restart !!! 191 192 ! in my changes here I assume that ht_0 AND r1_ht_0 have been updated 193 ! Note that the former calculation were using ht_0 so if it as not been updated ===>>> BUG 194 !!gm 195 196 188 197 ! Horizontal scale factor interpolations 189 198 ! -------------------------------------- 190 DO jk = 1,jpk 191 DO jj=1,jpj 192 DO ji=1,jpi 193 IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN 194 e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) / & 195 & ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) ) 196 ENDIF 197 END DO 198 END DO 199 END DO 200 ! 201 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 202 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 203 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 204 205 ! Vertical scale factor interpolations 206 ! ------------------------------------ 207 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 208 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 209 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 199 DO jj = 1, jpj 200 DO ji = 1, jpi 201 IF ( tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp ) THEN 202 DO jk = 1, jpk 203 e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) ) 204 END DO 205 ENDIF 206 END DO 207 END DO 208 ! 209 !!gm Note that if this routine is called in dom_vvl_init then all the lines below are uselss !!! 210 !! they are a duplication of dom_vvl_init lines 211 212 ! !== now fields ==! 213 ! 214 ! !* ssh at u- and v-points) 215 DO jj = 1, jpjm1 ; DO ji = 1, jpim1 ! start from 1 due to f-point 216 zsshu(ji,jj) = 0.5_wp * ( sshn(ji ,jj) + sshn(ji+1,jj ) ) * ssumask(ji,jj) 217 zsshv(ji,jj) = 0.5_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) ) * ssvmask(ji,jj) 218 zsshf(ji,jj) = 0.25_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 219 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) 220 END DO ; END DO 221 CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp , zsshf(:,:),'F', 1._wp ) 222 ! 223 ! !* hu and hv (and their inverse) 224 ht_n (:,:) = ht_0(:,:) + sshn(:,:) 225 hu_n (:,:) = hu_0(:,:) + zsshu(:,:) 226 hv_n (:,:) = hv_0(:,:) + zsshv(:,:) 227 r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) ! ss mask mask due to ISF 228 r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 229 ! 230 ! !* e3u, e3uw and e3v, e3vw 231 z_ssh_h0(:,:) = sshn (:,:) * r1_ht_0(:,:) ! t-point 232 DO jk = 1, jpkm1 233 e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * tmask(:,:,jk) ) 234 END DO 235 z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:) ! u-point 236 DO jk = 1, jpkm1 237 e3u_n (:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * umask(:,:,jk) ) 238 e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * wumask(:,:,jk) ) 239 END DO 240 z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:) ! v-point 241 DO jk = 1, jpkm1 242 e3v_n (:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * vmask(:,:,jk) ) 243 e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * wvmask(:,:,jk) ) 244 END DO 245 z_ssh_h0(:,:) = zsshf(:,:) * r1_hf_0(:,:) ! f-point 246 DO jk = 1, jpkm1 247 e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * fmask(:,:,jk) ) 248 END DO 249 250 z_ssh_h0(:,:) = 1._wp + sshn(:,:) * r1_ht_0(:,:) ! t-point 251 ! 252 IF( ln_isfcav ) THEN ! iceshelf cavities : ssh scaling not applied over the iceshelf thickness 253 DO jk = 1, jpkm1 254 gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:) 255 gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:) 256 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 257 END DO 258 ELSE 259 DO jk = 1, jpkm1 260 gdept_n(:,:,jk) = gdept_0(:,:,jk) * z_ssh_h0(:,:) 261 gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * z_ssh_h0(:,:) 262 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 263 END DO 264 ENDIF 210 265 211 ! t- and w- points depth 212 ! ---------------------- 213 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 214 gdepw_n(:,:,1) = 0.0_wp 215 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 216 DO jj = 1,jpj 217 DO ji = 1,jpi 218 DO jk = 2,mikt(ji,jj)-1 219 gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 220 gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 221 gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 222 END DO 223 IF (mikt(ji,jj) > 1) THEN 224 jk = mikt(ji,jj) 225 gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk) 226 gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 227 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk ) - sshn (ji,jj) 228 END IF 229 DO jk = mikt(ji,jj)+1, jpk 230 gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) 231 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 232 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk ) - sshn (ji,jj) 233 END DO 234 END DO 235 END DO 236 237 ! t-, u- and v- water column thickness 238 ! ------------------------------------ 239 ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp 240 DO jk = 1, jpkm1 241 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 242 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 243 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 244 END DO 245 ! ! Inverse of the local depth 246 r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 247 r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 248 249 END IF 266 ENDIF 250 267 251 268 !============================================================================= 252 269 ! compute velocity 253 270 ! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor). 254 ub(:,:,:) =un(:,:,:)255 vb(:,:,:) =vn(:,:,:)256 DO jk = 1, jpk257 DO jj = 1, jpj258 DO ji = 1, jpi271 ub(:,:,:) = un(:,:,:) 272 vb(:,:,:) = vn(:,:,:) 273 DO jk = 1, jpk 274 DO jj = 1, jpj 275 DO ji = 1, jpi 259 276 un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk); 260 277 vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk); … … 265 282 ! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column) 266 283 ! compute barotropic velocity now and after 267 ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 268 zbub(:,:) = SUM(ztrp,DIM=3) 269 ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 270 zbvb(:,:) = SUM(ztrp,DIM=3) 271 ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:); 272 zbun(:,:) = SUM(ztrp,DIM=3) 273 ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 274 zbvn(:,:) = SUM(ztrp,DIM=3) 284 zbub(:,:) = SUM( ub(:,:,:)*pe3u_b(:,:,:) , DIM=3 ) 285 zbvb(:,:) = SUM( vb(:,:,:)*pe3v_b(:,:,:) , DIM=3 ) 286 zbun(:,:) = SUM( un(:,:,:)* e3u_n(:,:,:) , DIM=3 ) 287 zbvn(:,:) = SUM( vn(:,:,:)* e3v_n(:,:,:) , DIM=3 ) 275 288 276 289 ! new water column … … 285 298 zucorr = 0._wp 286 299 zvcorr = 0._wp 287 DO jj = 1, jpj288 DO ji = 1, jpi300 DO jj = 1, jpj 301 DO ji = 1, jpi 289 302 IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN 290 303 zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj) … … 297 310 298 311 ! update velocity 299 DO jk = 1, jpk300 un(:,:,jk) =(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)301 vn(:,:,jk) =(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)312 DO jk = 1, jpk 313 un(:,:,jk) = ( un(:,:,jk) - zucorr(:,:) ) * umask(:,:,jk) 314 vn(:,:,jk) = ( vn(:,:,jk) - zvcorr(:,:) ) * vmask(:,:,jk) 302 315 END DO 303 316 … … 305 318 ! compute temp and salt 306 319 ! compute new tn and sn if we open a new cell 307 tsb (:,:,:,:) = tsn(:,:,:,:)308 zts0 (:,:,:,:) = tsn(:,:,:,:)320 tsb (:,:,:,:) = tsn (:,:,:,:) 321 zts0 (:,:,:,:) = tsn (:,:,:,:) 309 322 ztmask1(:,:,:) = ptmask_b(:,:,:) 310 323 ztmask0(:,:,:) = ptmask_b(:,:,:) 311 324 DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case) 312 313 314 315 316 317 318 319 320 321 322 & 323 & 324 & 325 326 & 327 & 328 & 329 330 331 332 333 334 335 336 & +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk337 338 & +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk339 340 ENDIF341 ENDIF342 343 344 345 346 347 348 349 350 351 325 DO jk = 1,jpk-1 326 zdmask=tmask(:,:,jk)-ztmask0(:,:,jk); 327 DO jj = 2,jpj-1 328 DO ji = fs_2,fs_jpim1 329 jip1=ji+1; jim1=ji-1; 330 jjp1=jj+1; jjm1=jj-1; 331 summsk= (ztmask0(jip1,jj ,jk)+ztmask0(jim1,jj ,jk)+ztmask0(ji ,jjp1,jk)+ztmask0(ji ,jjm1,jk)) 332 IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN 333 !! horizontal basic extrapolation 334 tsn(ji,jj,jk,1)=( zts0(jip1,jj ,jk,1)*ztmask0(jip1,jj ,jk) & 335 & +zts0(jim1,jj ,jk,1)*ztmask0(jim1,jj ,jk) & 336 & +zts0(ji ,jjp1,jk,1)*ztmask0(ji ,jjp1,jk) & 337 & +zts0(ji ,jjm1,jk,1)*ztmask0(ji ,jjm1,jk) ) / summsk 338 tsn(ji,jj,jk,2)=( zts0(jip1,jj ,jk,2)*ztmask0(jip1,jj ,jk) & 339 & +zts0(jim1,jj ,jk,2)*ztmask0(jim1,jj ,jk) & 340 & +zts0(ji ,jjp1,jk,2)*ztmask0(ji ,jjp1,jk) & 341 & +zts0(ji ,jjm1,jk,2)*ztmask0(ji ,jjm1,jk) ) / summsk 342 ztmask1(ji,jj,jk)=1 343 ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN 344 !! vertical extrapolation if horizontal extrapolation failed 345 jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1) 346 summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1)) 347 IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN 348 tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1) & 349 & +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / summsk 350 tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1) & 351 & +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / summsk 352 ztmask1(ji,jj,jk)=1._wp 353 ENDIF 354 ENDIF 355 END DO 356 END DO 357 END DO 358 ! 359 CALL lbc_lnk_multi( tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.) 360 ! 361 ! update 362 zts0(:,:,:,:) = tsn(:,:,:,:) 363 ztmask0 = ztmask1 364 ! 352 365 END DO 353 366 … … 395 408 ! set mbkt and mikt to 1 in thiese location 396 409 WHERE (SUM(tmask,dim=3) == 0) 397 mbkt(:,:) =1 ; mbku(:,:)=1 ; mbkv(:,:)=1398 mikt(:,:) =1 ; miku(:,:)=1 ; mikv(:,:)=1410 mbkt(:,:) = 1 ; mbku(:,:)=1 ; mbkv(:,:) = 1 411 mikt(:,:) = 1 ; miku(:,:)=1 ; mikv(:,:) = 1 399 412 END WHERE 400 413 ! -------------------------------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.