Changeset 9987 for branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
- Timestamp:
- 2018-07-23T11:33:03+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r7960 r9987 26 26 USE zdfbfr 27 27 USE fldread ! read input field at current time step 28 29 28 USE lib_fortran, ONLY: glob_sum 30 29 31 30 IMPLICIT NONE … … 53 52 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: risfLeff !:effective length (Leff) BG03 nn_isf==2 54 53 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point 55 #if defined key_agrif56 ! AGRIF can not handle these arrays as integers. The reason is a mystery but problems avoided by declaring them as reals57 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base58 !: (first wet level and last level include in the tbl)59 #else60 54 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base 61 #endif62 55 63 56 … … 90 83 INTEGER :: ji, jj, jk, ijkmin, inum, ierror 91 84 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 85 REAL(wp) :: zgreenland_fwfisf_sum, zantarctica_fwfisf_sum 92 86 REAL(wp) :: rmin 93 87 REAL(wp) :: zhk 94 CHARACTER(len=256) :: cfisf, cvarzisf, cvarhisf ! name for isf file 88 REAL(wp) :: zt_frz, zpress 89 CHARACTER(len=256) :: cfisf , cvarzisf, cvarhisf ! name for isf file 95 90 CHARACTER(LEN=256) :: cnameis ! name of iceshelf file 96 91 CHARACTER (LEN=32) :: cvarLeff ! variable name for efficient Length scale 97 92 INTEGER :: ios ! Local integer output status for namelist read 93 94 REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d 95 REAL(wp), DIMENSION(:,: ), POINTER :: zqhcisf2d 98 96 ! 99 97 !!--------------------------------------------------------------------- … … 176 174 DO jj = 1, jpj 177 175 jk = 2 178 DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ; jk = jk + 1 ; END DO176 DO WHILE ( jk .LE. mbkt(ji,jj) .AND. gdepw_0(ji,jj,jk) < rzisf_tbl(ji,jj) ) ; jk = jk + 1 ; END DO 179 177 misfkt(ji,jj) = jk-1 180 178 END DO … … 194 192 END IF 195 193 194 ! save initial top boundary layer thickness 196 195 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 196 197 END IF 198 199 ! ! ---------------------------------------- ! 200 IF( kt /= nit000 ) THEN ! Swap of forcing fields ! 201 ! ! ---------------------------------------- ! 202 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000 203 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine 204 ! 205 ENDIF 206 207 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 197 208 198 209 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf … … 205 216 206 217 ! determine the deepest level influenced by the boundary layer 207 ! test on tmask useless ?????208 218 DO jk = ikt, mbkt(ji,jj) 209 219 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk … … 217 227 END DO 218 228 END DO 219 220 END IF221 222 ! ! ---------------------------------------- !223 IF( kt /= nit000 ) THEN ! Swap of forcing fields !224 ! ! ---------------------------------------- !225 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000226 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine227 !228 ENDIF229 230 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN231 232 229 233 230 ! compute salf and heat flux … … 256 253 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 257 254 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) 255 256 IF( lk_oasis) THEN 257 ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true 258 IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN 259 260 ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern 261 ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets 262 ! to preserve total freshwater conservation in coupled models without an active ice sheet model. 263 264 ! All related global sums must be done bit reproducibly 265 zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 266 267 ! use ABS function because we need to preserve the sign of fwfisf 268 WHERE( greenland_icesheet_mask(:,:) == 1.0 ) & 269 & fwfisf(:,:) = fwfisf(:,:) * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) & 270 & / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) ) 271 272 ! check 273 IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum 274 275 zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 276 277 IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum 278 279 zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 280 281 ! use ABS function because we need to preserve the sign of fwfisf 282 WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) & 283 & fwfisf(:,:) = fwfisf(:,:) * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) & 284 & / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) ) 285 286 ! check 287 IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum 288 289 zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 290 291 IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum 292 293 ENDIF 294 ENDIF 295 258 296 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 259 297 stbl(:,:) = soce … … 264 302 !CALL fld_read ( kt, nn_fsbc, sf_qisf ) 265 303 fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1) ! fwf 304 305 IF( lk_oasis) THEN 306 ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true 307 IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN 308 309 ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern 310 ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets 311 ! to preserve total freshwater conservation in coupled models without an active ice sheet model. 312 313 ! All related global sums must be done bit reproducibly 314 zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 315 316 ! use ABS function because we need to preserve the sign of fwfisf 317 WHERE( greenland_icesheet_mask(:,:) == 1.0 ) & 318 & fwfisf(:,:) = fwfisf(:,:) * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) & 319 & / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) ) 320 321 ! check 322 IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum 323 324 zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 325 326 IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum 327 328 zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 329 330 ! use ABS function because we need to preserve the sign of fwfisf 331 WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) & 332 & fwfisf(:,:) = fwfisf(:,:) * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) & 333 & / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) ) 334 335 ! check 336 IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum 337 338 zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 339 340 IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum 341 342 ENDIF 343 ENDIF 344 266 345 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 267 346 !qisf(:,:) = sf_qisf(1)%fnow(:,:,1) ! heat flux … … 270 349 END IF 271 350 ! compute tsc due to isf 272 ! WARNING water add at temp = 0C, correction term is added in trasbc, maybe better here but need a 3D variable). 273 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp ! 351 ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable). 352 ! zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 353 zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 354 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 ! 274 355 275 356 ! salt effect already take into account in vertical advection 276 357 risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0 277 358 359 ! output 360 IF( iom_use('qlatisf' ) ) CALL iom_put('qlatisf', qisf) 361 IF( iom_use('fwfisf' ) ) CALL iom_put('fwfisf' , fwfisf * stbl(:,:) / soce ) 362 363 ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now 364 fwfisf(:,:) = rdivisf * fwfisf(:,:) 365 278 366 ! lbclnk 279 367 CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) … … 281 369 CALL lbc_lnk(fwfisf(:,:) ,'T',1.) 282 370 CALL lbc_lnk(qisf(:,:) ,'T',1.) 371 372 !============================================================================================================================================= 373 IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 374 CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 375 CALL wrk_alloc( jpi,jpj, zqhcisf2d ) 376 377 zfwfisf3d(:,:,:) = 0.0_wp ! 3d ice shelf melting (kg/m2/s) 378 zqhcisf3d(:,:,:) = 0.0_wp ! 3d heat content flux (W/m2) 379 zqlatisf3d(:,:,:)= 0.0_wp ! 3d ice shelf melting latent heat flux (W/m2) 380 zqhcisf2d(:,:) = fwfisf(:,:) * zt_frz * rcp ! 2d heat content flux (W/m2) 381 382 DO jj = 1,jpj 383 DO ji = 1,jpi 384 ikt = misfkt(ji,jj) 385 ikb = misfkb(ji,jj) 386 DO jk = ikt, ikb - 1 387 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 388 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 389 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 390 END DO 391 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 392 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 393 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 394 END DO 395 END DO 396 397 CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:)) 398 CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:)) 399 CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:)) 400 CALL iom_put('qhcisf' , zqhcisf2d (:,: )) 401 402 CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 403 CALL wrk_dealloc( jpi,jpj, zqhcisf2d ) 404 END IF 405 !============================================================================================================================================= 283 406 284 407 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! … … 295 418 ENDIF 296 419 ! 297 ! output298 CALL iom_put('qisf' , qisf)299 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )300 420 END IF 301 421 … … 370 490 ! Calculate freezing temperature 371 491 zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 372 zt_frz = eos_fzp(tsb(ji,jj,ik,jp_sal), zpress)492 CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress) 373 493 zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik) ! sum temp 374 494 ENDDO … … 452 572 zti(:,:)=tinsitu( ttbl, stbl, zpress ) 453 573 ! Calculate freezing temperature 454 zfrz(:,:)=eos_fzp( sss_m(:,:), zpress )574 CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress ) 455 575 456 576 … … 472 592 473 593 nit = nit + 1 474 IF (nit .GE. 100) THEN 475 !WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b,zgammat, rn_gammat0, rn_tfri2, nn_gammablk, ji,jj 476 !WRITE(numout,*) "sbcisf : too many iteration ... ", (zhtflx - zhtflx_b)/zhtflx 477 CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 478 END IF 594 IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 595 479 596 ! save gammat and compute zhtflx_b 480 597 zgammat2d(ji,jj)=zgammat … … 794 911 ! test on tmask useless ????? 795 912 DO jk = ikt, mbkt(ji,jj) 796 !IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk913 IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 797 914 END DO 798 915 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
Note: See TracChangeset
for help on using the changeset viewer.
