- Timestamp:
- 2017-03-17T08:46:30+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r7256 r7806 32 32 PRIVATE 33 33 34 PUBLIC sbc_isf, sbc_isf_ div, sbc_isf_alloc ! routine called in sbcmod and divcur34 PUBLIC sbc_isf, sbc_isf_init, sbc_isf_div, sbc_isf_alloc ! routine called in sbcmod and divcur 35 35 36 36 ! public in order to be able to output then … … 54 54 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point 55 55 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base 56 57 LOGICAL, PUBLIC :: l_isfcpl = .false. ! isf recieved from oasis 56 58 57 59 … … 81 83 82 84 SUBROUTINE sbc_isf(kt) 85 83 86 INTEGER, INTENT(in) :: kt ! ocean time step 87 INTEGER :: ji, jj, jk 88 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 89 REAL(wp) :: zhk 90 REAL(wp) :: zt_frz, zpress 91 REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d 92 REAL(wp), DIMENSION(:,: ), POINTER :: zqhcisf2d 93 REAL(wp) :: zhisf 94 95 96 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 97 98 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 99 DO jj = 1,jpj 100 DO ji = 1,jpi 101 ikt = misfkt(ji,jj) 102 ikb = misfkt(ji,jj) 103 ! thickness of boundary layer at least the top level thickness 104 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt)) 105 106 ! determine the deepest level influenced by the boundary layer 107 DO jk = ikt, mbkt(ji,jj) 108 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 109 END DO 110 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 111 misfkb(ji,jj) = ikb ! last wet level of the tbl 112 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 113 114 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 115 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 116 END DO 117 END DO 118 119 ! compute salf and heat flux 120 IF (nn_isf == 1) THEN 121 ! realistic ice shelf formulation 122 ! compute T/S/U/V for the top boundary layer 123 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 124 CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 125 CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U') 126 CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V') 127 ! iom print 128 CALL iom_put('ttbl',ttbl(:,:)) 129 CALL iom_put('stbl',stbl(:,:)) 130 CALL iom_put('utbl',utbl(:,:)) 131 CALL iom_put('vtbl',vtbl(:,:)) 132 ! compute fwf and heat flux 133 IF( .NOT.l_isfcpl ) THEN ; CALL sbc_isf_cav (kt) 134 ELSE ; qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 135 ENDIF 136 137 ELSE IF (nn_isf == 2) THEN 138 ! Beckmann and Goosse parametrisation 139 stbl(:,:) = soce 140 CALL sbc_isf_bg03(kt) 141 142 ELSE IF (nn_isf == 3) THEN 143 ! specified runoff in depth (Mathiot et al., XXXX in preparation) 144 IF( .NOT.l_isfcpl ) THEN 145 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 146 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) 147 ENDIF 148 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 149 stbl(:,:) = soce 150 151 ELSE IF (nn_isf == 4) THEN 152 ! specified fwf and heat flux forcing beneath the ice shelf 153 IF( .NOT.l_isfcpl ) THEN 154 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 155 !CALL fld_read ( kt, nn_fsbc, sf_qisf ) 156 fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1) ! fwf 157 ENDIF 158 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 159 !qisf(:,:) = sf_qisf(1)%fnow(:,:,1) ! heat flux 160 stbl(:,:) = soce 161 162 END IF 163 ! compute tsc due to isf 164 ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable). 165 ! zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 166 zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 167 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 ! 168 169 ! salt effect already take into account in vertical advection 170 risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0 171 172 ! output 173 IF( iom_use('qlatisf' ) ) CALL iom_put('qlatisf', qisf) 174 IF( iom_use('fwfisf' ) ) CALL iom_put('fwfisf' , fwfisf * stbl(:,:) / soce ) 175 176 ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now 177 fwfisf(:,:) = rdivisf * fwfisf(:,:) 178 179 ! lbclnk 180 CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) 181 CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.) 182 CALL lbc_lnk(fwfisf(:,:) ,'T',1.) 183 CALL lbc_lnk(qisf(:,:) ,'T',1.) 184 185 ! Diagnostics 186 IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 187 ! 188 CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 189 CALL wrk_alloc( jpi,jpj, zqhcisf2d ) 190 ! 191 zfwfisf3d(:,:,:) = 0.0_wp ! 3d ice shelf melting (kg/m2/s) 192 zqhcisf3d(:,:,:) = 0.0_wp ! 3d heat content flux (W/m2) 193 zqlatisf3d(:,:,:)= 0.0_wp ! 3d ice shelf melting latent heat flux (W/m2) 194 zqhcisf2d(:,:) = fwfisf(:,:) * zt_frz * rcp ! 2d heat content flux (W/m2) 195 ! 196 DO jj = 1,jpj 197 DO ji = 1,jpi 198 ikt = misfkt(ji,jj) 199 ikb = misfkb(ji,jj) 200 DO jk = ikt, ikb - 1 201 zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 202 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf(ji,jj) * zhisf 203 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf 204 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf(ji,jj) * zhisf 205 END DO 206 jk = ikb 207 zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 208 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * zhisf * ralpha(ji,jj) 209 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf * ralpha(ji,jj) 210 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * zhisf * ralpha(ji,jj) 211 END DO 212 END DO 213 ! 214 CALL iom_put( 'fwfisf3d' , zfwfisf3d (:,:,:) ) 215 CALL iom_put( 'qlatisf3d', zqlatisf3d(:,:,:) ) 216 CALL iom_put( 'qhcisf3d' , zqhcisf3d (:,:,:) ) 217 CALL iom_put( 'qhcisf' , zqhcisf2d (:,: ) ) 218 ! 219 CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 220 CALL wrk_dealloc( jpi,jpj, zqhcisf2d ) 221 ! 222 END IF 223 ! 224 END IF 225 ! 226 ! 227 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! 228 IF( ln_rstart .AND. & ! Restart: read in restart file 229 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN 230 IF(lwp) WRITE(numout,*) ' nit000-1 isf tracer content forcing fields read in the restart file' 231 CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) ) ! before salt content isf_tsc trend 232 CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) ) ! before salt content isf_tsc trend 233 CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) ) ! before salt content isf_tsc trend 234 ELSE 235 fwfisf_b(:,:) = fwfisf(:,:) 236 risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 237 END IF 238 ENDIF 239 ! 240 IF( lrst_oce ) THEN 241 IF(lwp) WRITE(numout,*) 242 IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ', & 243 & 'at it= ', kt,' date= ', ndastp 244 IF(lwp) WRITE(numout,*) '~~~~' 245 CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) ) 246 CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) ) 247 CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) ) 248 ENDIF 249 ! 250 END SUBROUTINE sbc_isf 251 252 SUBROUTINE sbc_isf_init 253 84 254 INTEGER :: ji, jj, jk, ijkmin, inum, ierror 85 255 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 86 REAL(wp) :: rmin87 256 REAL(wp) :: zhk 88 REAL(wp) :: zt_frz, zpress89 257 CHARACTER(len=256) :: cfisf , cvarzisf, cvarhisf ! name for isf file 90 258 CHARACTER(LEN=256) :: cnameis ! name of iceshelf file 91 259 CHARACTER (LEN=32) :: cvarLeff ! variable name for efficient Length scale 92 260 INTEGER :: ios ! Local integer output status for namelist read 261 93 262 ! 94 263 !!--------------------------------------------------------------------- … … 97 266 ! 98 267 ! 99 ! ! ====================== !100 IF( kt == nit000 ) THEN ! First call kt=nit000 !101 ! ! ====================== !102 268 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 103 269 READ ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) … … 139 305 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 140 306 ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN 141 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 142 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 143 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 307 IF( .NOT.l_isfcpl ) THEN 308 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 309 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 310 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 311 ENDIF 144 312 145 313 !: read effective lenght (BG03) … … 182 350 183 351 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 184 ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 185 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 186 ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) ) 187 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 188 !CALL fld_fill( sf_qisf , (/ sn_qisf /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data' , 'namsbc_isf' ) 352 IF( .NOT.l_isfcpl ) THEN 353 ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 354 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 355 ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) ) 356 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 357 !CALL fld_fill( sf_qisf , (/ sn_qisf /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data' , 'namsbc_isf' ) 358 ENDIF 189 359 END IF 190 191 360 ! save initial top boundary layer thickness 192 361 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 193 194 END IF195 196 ! ! ---------------------------------------- !197 IF( kt /= nit000 ) THEN ! Swap of forcing fields !198 ! ! ---------------------------------------- !199 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000200 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine201 !202 ENDIF203 204 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN205 206 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf207 DO jj = 1,jpj208 DO ji = 1,jpi209 ikt = misfkt(ji,jj)210 ikb = misfkt(ji,jj)211 ! thickness of boundary layer at least the top level thickness212 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))213 214 ! determine the deepest level influenced by the boundary layer215 DO jk = ikt, mbkt(ji,jj)216 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk217 END DO218 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.219 misfkb(ji,jj) = ikb ! last wet level of the tbl220 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)221 222 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1223 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer224 END DO225 END DO226 227 ! compute salf and heat flux228 IF (nn_isf == 1) THEN229 ! realistic ice shelf formulation230 ! compute T/S/U/V for the top boundary layer231 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')232 CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')233 CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U')234 CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V')235 ! iom print236 CALL iom_put('ttbl',ttbl(:,:))237 CALL iom_put('stbl',stbl(:,:))238 CALL iom_put('utbl',utbl(:,:))239 CALL iom_put('vtbl',vtbl(:,:))240 ! compute fwf and heat flux241 CALL sbc_isf_cav (kt)242 243 ELSE IF (nn_isf == 2) THEN244 ! Beckmann and Goosse parametrisation245 stbl(:,:) = soce246 CALL sbc_isf_bg03(kt)247 248 ELSE IF (nn_isf == 3) THEN249 ! specified runoff in depth (Mathiot et al., XXXX in preparation)250 CALL fld_read ( kt, nn_fsbc, sf_rnfisf )251 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting)252 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux253 stbl(:,:) = soce254 255 ELSE IF (nn_isf == 4) THEN256 ! specified fwf and heat flux forcing beneath the ice shelf257 CALL fld_read ( kt, nn_fsbc, sf_fwfisf )258 !CALL fld_read ( kt, nn_fsbc, sf_qisf )259 fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1) ! fwf260 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux261 !qisf(:,:) = sf_qisf(1)%fnow(:,:,1) ! heat flux262 stbl(:,:) = soce263 264 END IF265 ! compute tsc due to isf266 ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable).267 ! zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04268 zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )269 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 !270 271 ! salt effect already take into account in vertical advection272 risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0273 274 ! output275 IF( iom_use('qisf' ) ) CALL iom_put('qisf' , qisf)276 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )277 278 ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now279 fwfisf(:,:) = rdivisf * fwfisf(:,:)280 281 ! lbclnk282 CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)283 CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)284 CALL lbc_lnk(fwfisf(:,:) ,'T',1.)285 CALL lbc_lnk(qisf(:,:) ,'T',1.)286 287 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 !288 IF( ln_rstart .AND. & ! Restart: read in restart file289 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN290 IF(lwp) WRITE(numout,*) ' nit000-1 isf tracer content forcing fields read in the restart file'291 CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) ) ! before salt content isf_tsc trend292 CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) ) ! before salt content isf_tsc trend293 CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) ) ! before salt content isf_tsc trend294 ELSE295 fwfisf_b(:,:) = fwfisf(:,:)296 risf_tsc_b(:,:,:)= risf_tsc(:,:,:)297 END IF298 ENDIF299 362 ! 300 END IF301 302 END SUBROUTINE sbc_isf 363 END SUBROUTINE sbc_isf_init 364 365 303 366 304 367 INTEGER FUNCTION sbc_isf_alloc()
Note: See TracChangeset
for help on using the changeset viewer.
