Changeset 11395
- Timestamp:
- 2019-08-02T16:19:00+02:00 (6 years ago)
- Location:
- NEMO/branches/2019/ENHANCE-02_ISF_nemo
- Files:
-
- 14 added
- 4 deleted
- 25 edited
- 5 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/ENHANCE-02_ISF_nemo/cfgs/SHARED/field_def_nemo-oce.xml
r10824 r11395 264 264 265 265 <!-- * variable related to ice shelf forcing * --> 266 <field id="fwfisf" long_name="Ice shelf melting" unit="kg/m2/s" /> 267 <field id="fwfisf3d" long_name="Ice shelf melting" unit="kg/m2/s" grid_ref="grid_T_3D" /> 268 <field id="qlatisf" long_name="Ice shelf latent heat flux" unit="W/m2" /> 269 <field id="qlatisf3d" long_name="Ice shelf latent heat flux" unit="W/m2" grid_ref="grid_T_3D" /> 270 <field id="qhcisf" long_name="Ice shelf heat content flux" unit="W/m2" /> 271 <field id="qhcisf3d" long_name="Ice shelf heat content flux" unit="W/m2" grid_ref="grid_T_3D" /> 272 <field id="isfgammat" long_name="transfert coefficient for isf (temperature) " unit="m/s" /> 273 <field id="isfgammas" long_name="transfert coefficient for isf (salinity) " unit="m/s" /> 274 <field id="stbl" long_name="salinity in the Losh tbl " unit="PSU" /> 275 <field id="ttbl" long_name="temperature in the Losh tbl " unit="C" /> 276 <field id="utbl" long_name="zonal current in the Losh tbl at T point " unit="m/s" /> 277 <field id="vtbl" long_name="merid current in the Losh tbl at T point " unit="m/s" /> 278 <field id="thermald" long_name="thermal driving of ice shelf melting " unit="C" /> 279 <field id="tfrz" long_name="top freezing point (used to compute melt) " unit="C" /> 280 <field id="tinsitu" long_name="top insitu temperature (used to cmpt melt) " unit="C" /> 281 <field id="ustar" long_name="ustar at T point used in ice shelf melting " unit="m/s" /> 266 <field id="fwfisf_cav" long_name="Ice shelf melting" unit="kg/m2/s" /> 267 <field id="fwfisf_par" long_name="Ice shelf melting" unit="kg/m2/s" /> 268 <field id="qoceisf_cav" long_name="Ice shelf ocean heat flux" unit="W/m2" /> 269 <field id="qoceisf_par" long_name="Ice shelf ocean heat flux" unit="W/m2" /> 270 <field id="qlatisf_cav" long_name="Ice shelf latent heat flux" unit="W/m2" /> 271 <field id="qlatisf_par" long_name="Ice shelf latent heat flux" unit="W/m2" /> 272 <field id="qhcisf_cav" long_name="Ice shelf heat content flux" unit="W/m2" /> 273 <field id="qhcisf_par" long_name="Ice shelf heat content flux" unit="W/m2" /> 274 <field id="fwfisf3d_cav" long_name="Ice shelf melting" unit="kg/m2/s" grid_ref="grid_T_3D" /> 275 <field id="fwfisf3d_par" long_name="Ice shelf melting" unit="kg/m2/s" grid_ref="grid_T_3D" /> 276 <field id="qoceisf3d_cav" long_name="Ice shelf ocean heat flux" unit="W/m2" grid_ref="grid_T_3D" /> 277 <field id="qoceisf3d_par" long_name="Ice shelf ocean heat flux" unit="W/m2" grid_ref="grid_T_3D" /> 278 <field id="qlatisf3d_cav" long_name="Ice shelf latent heat flux" unit="W/m2" grid_ref="grid_T_3D" /> 279 <field id="qlatisf3d_par" long_name="Ice shelf latent heat flux" unit="W/m2" grid_ref="grid_T_3D" /> 280 <field id="qhcisf3d_cav" long_name="Ice shelf heat content flux" unit="W/m2" grid_ref="grid_T_3D" /> 281 <field id="qhcisf3d_par" long_name="Ice shelf heat content flux" unit="W/m2" grid_ref="grid_T_3D" /> 282 <field id="ttbl_cav" long_name="temperature in the tracer sample depth " unit="C" /> 283 <field id="ttbl_par" long_name="temperature in the tracer sample depth " unit="C" /> 284 <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting " unit="C" /> 285 <field id="isfthermald_par" long_name="thermal driving of ice shelf melting " unit="C" /> 286 <field id="isfgammat" long_name="transfert coefficient for isf (temperature) " unit="m/s" /> 287 <field id="isfgammas" long_name="transfert coefficient for isf (salinity) " unit="m/s" /> 288 <field id="stbl" long_name="salinity in the Losh tbl " unit="PSU" /> 289 <field id="utbl" long_name="zonal current in the Losh tbl at T point " unit="m/s" /> 290 <field id="vtbl" long_name="merid current in the Losh tbl at T point " unit="m/s" /> 291 <field id="isfustar" long_name="ustar at T point used in ice shelf melting " unit="m/s" /> 292 <field id="qconisf" long_name="Conductive heat flux through the ice shelf" unit="W/m2" /> 282 293 283 294 <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/cfgs/SHARED/namelist_ref
r10808 r11395 52 52 cn_ocerst_outdir = "." ! directory in which to write output ocean restarts 53 53 ln_iscpl = .false. ! cavity evolution forcing or coupling to ice sheet model 54 nn_istate = 0! output the initial state (1) or not (0)54 nn_istate = 1 ! output the initial state (1) or not (0) 55 55 ln_rst_list = .false. ! output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 56 56 nn_stock = 5840 ! frequency of creation of a restart file (modulo referenced to 1) … … 68 68 !----------------------------------------------------------------------- 69 69 ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time 70 rn_isfhmin = 1.00 ! treshold [m] to discriminate grounding ice from floating ice71 70 ! 72 71 rn_rdt = 5400. ! time step for the dynamics and tracer … … 75 74 ln_crs = .false. ! Logical switch for coarsening module (T => fill namcrs) 76 75 ! 77 ln_meshmask = . false. ! =T create a mesh file76 ln_meshmask = .true. ! =T create a mesh file 78 77 / 79 78 !----------------------------------------------------------------------- … … 180 179 !! namsbc_rnf river runoffs (ln_rnf =T) 181 180 !! namsbc_apr Atmospheric Pressure (ln_apr_dyn =T) 182 !! namsbc_isf ice shelf melting/freezing (ln_isfcav =T : read (ln_read_cfg=T) or set or usr_def_zgr )183 181 !! namsbc_iscpl coupling option between land ice model and ocean (ln_isfcav =T) 184 182 !! namsbc_wave external fields from wave model (ln_wave =T) … … 436 434 / 437 435 !----------------------------------------------------------------------- 438 &nam sbc_isf ! Top boundary layer (ISF) (ln_isfcav =T: read (ln_read_cfg=T)439 !----------------------------------------------------------------------- or set or usr_def_zgr )436 &namisf ! Top boundary layer (ISF) (ln_isf = T .AND. ln_isfcav: read (ln_read_cfg=T) 437 !----------------------------------------------------------------------- or set or usr_def_zgr ) 440 438 ! ! type of top boundary layer 441 nn_isf = 1 ! ice shelf melting/freezing 442 ! 1 = presence of ISF ; 2 = bg03 parametrisation 443 ! 3 = rnf file for ISF ; 4 = ISF specified freshwater flux 444 ! options 1 and 4 need ln_isfcav = .true. (domzgr) 445 ! ! nn_isf = 1 or 2 cases: 439 ln_isfcav_mlt = .false. ! ice shelf melting into the cavity 440 cn_isfcav_mlt = '3eq' ! ice shelf melting formulation (spe/2eq/3eq/oasis) 441 ! ! spe = fwfisf is read from a forcing field 442 ! ! 2eq = ISOMIP like: 2 equations formulation (Hunter et al., 2006) 443 ! ! 3eq = ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015) 444 ! ! oasis = fwfisf is given by oasis 445 ! ! cn_isfcav_mlt = 2eq or 3eq cases: 446 cn_gammablk = 'ad15' ! scheme to compute gammat/s (spe,ad15,hj99) 447 ! ! ad15 = velocity dependend Gamma (u* * gammat/s) (Jenkins et al. 2010) 448 ! ! hj99 = velocity and stability dependent Gamma (Holland et al. 1999) 446 449 rn_gammat0 = 1.e-4 ! gammat coefficient used in blk formula 447 450 rn_gammas0 = 1.e-4 ! gammas coefficient used in blk formula 448 ! ! nn_isf = 1 or 4 cases:449 rn_h isf_tbl= 30. ! thickness of the top boundary layer (Losh et al. 2008)451 ! 452 rn_htbl = 30. ! thickness of the top boundary layer (Losh et al. 2008) 450 453 ! ! 0 => thickness of the tbl = thickness of the first wet cell 451 ! ! nn_isf = 1 case 452 nn_isfblk = 1 ! 1 ISOMIP like: 2 equations formulation (Hunter et al., 2006) 453 ! ! 2 ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015) 454 nn_gammablk = 1 ! 0 = cst Gammat (= gammat/s) 455 ! ! 1 = velocity dependend Gamma (u* * gammat/s) (Jenkins et al. 2010) 456 ! ! 2 = velocity and stability dependent Gamma (Holland et al. 1999) 457 458 !___________!_____________!___________________!___________!_____________!_________!___________!__________!__________!_______________! 459 ! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 460 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 461 !* nn_isf = 4 case 462 sn_fwfisf = 'rnfisf' , -12 ,'sowflisf' , .false. , .true. , 'yearly' , '' , '' , '' 463 !* nn_isf = 3 case 464 sn_rnfisf = 'rnfisf' , -12 ,'sofwfisf' , .false. , .true. , 'yearly' , '' , '' , '' 465 !* nn_isf = 2 and 3 cases 466 sn_depmax_isf ='rnfisf' , -12 ,'sozisfmax', .false. , .true. , 'yearly' , '' , '' , '' 467 sn_depmin_isf ='rnfisf' , -12 ,'sozisfmin', .false. , .true. , 'yearly' , '' , '' , '' 468 !* nn_isf = 2 case 469 sn_Leff_isf = 'rnfisf' , -12 ,'Leff' , .false. , .true. , 'yearly' , '' , '' , '' 454 ! 455 !* 'spe' and 'oasis' case 456 !___________!_____________!___________________!___________!_____________!_________!___________!__________!__________!_______________! 457 ! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 458 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 459 sn_isfcav_fwf = 'isfmlt_cav', -12 , 'fwflisf' , .false. , .true. , 'yearly' , '' , '' , '' 460 ! 461 ! 462 ln_isfpar_mlt = .false. ! ice shelf melting parametrised 463 cn_isfpar_mlt = 'spe' ! ice shelf melting parametrisation (spe/bg03/oasis) 464 ! 465 !* all cases 466 !___________!_____________!___________________!___________!_____________!_________!___________!__________!__________!_______________! 467 ! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 468 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 469 sn_isfpar_zmax = 'isfmlt_par', 0 ,'sozisfmax', .false. , .true. , 'yearly' , '' , '' , '' 470 sn_isfpar_zmin = 'isfmlt_par', 0 ,'sozisfmin', .false. , .true. , 'yearly' , '' , '' , '' 471 !* 'spe' and 'oasis' case 472 sn_isfpar_fwf = 'isfmlt_par' , -12 ,'fwfisf' , .false. , .true. , 'yearly' , '' , '' , '' 473 !* 'bg03' case 474 sn_isfpar_Leff = 'isfmlt_par', 0 ,'Leff' , .false. , .true. , 'yearly' , '' , '' , '' 475 ! 476 ! 477 !ln_iscpl = .false. 478 ! nn_drown = 10 ! number of iteration of the extrapolation loop (fill the new wet cells) 479 ! ln_hsb = .false. ! activate conservation module (conservation exact after a time of rn_fiscpl) 480 ! nn_fiscpl = 43800 ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency) 470 481 / 471 482 !----------------------------------------------------------------------- -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/BDY/bdyvol.F90
r10481 r11395 16 16 USE dom_oce ! ocean space and time domain 17 17 USE phycst ! physical constants 18 USE sbcisf! ice shelf18 USE isf ! ice shelf 19 19 ! 20 20 USE in_out_manager ! I/O manager … … 77 77 ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 78 78 ! ----------------------------------------------------------------------- 79 IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf (:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau079 IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0 80 80 81 81 ! Compute bdy surface each cycle if non linear free surface -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA/diahsb.F90
r10425 r11395 18 18 USE sbc_oce ! surface thermohaline fluxes 19 19 USE sbcrnf ! river runoff 20 USE sbcisf! ice shelves20 USE isf ! ice shelves 21 21 USE domvvl ! vertical scale factors 22 22 USE traqsr ! penetrative solar radiation … … 91 91 ! 1 - Trends due to forcing ! 92 92 ! ------------------------- ! 93 z_frc_trd_v = r1_rau0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf (:,:) ) * surf(:,:) ) ! volume fluxes93 z_frc_trd_v = r1_rau0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) ) ! volume fluxes 94 94 z_frc_trd_t = glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) ) ! heat fluxes 95 95 z_frc_trd_s = glob_sum( 'diahsb', sbc_tsc(:,:,jp_sal) * surf(:,:) ) ! salt fluxes … … 98 98 IF( ln_rnf_sal) z_frc_trd_s = z_frc_trd_s + glob_sum( 'diahsb', rnf_tsc(:,:,jp_sal) * surf(:,:) ) 99 99 ! ! Add ice shelf heat & salt input 100 IF( ln_isf ) z_frc_trd_t = z_frc_trd_t + glob_sum( 'diahsb', risf_tsc(:,:,jp_tem) * surf(:,:) ) 100 IF( ln_isf ) z_frc_trd_t = z_frc_trd_t & 101 & + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) ) 101 102 ! ! Add penetrative solar radiation 102 103 IF( ln_traqsr ) z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( 'diahsb', qsr (:,:) * surf(:,:) ) -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90
r10425 r11395 26 26 !!---------------------------------------------------------------------- 27 27 USE oce ! ocean dynamics and tracers 28 USE isf 28 29 USE dom_oce ! ocean space and time domain 29 30 USE phycst ! physical constants … … 904 905 CALL iom_rstput( 0, 0, inum, 'vomecrty', vn ) ! now j-velocity 905 906 CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity 907 CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep ) ! now k-velocity 908 CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav ) ! now k-velocity 909 CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav ) ! now k-velocity 910 CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav ) ! now k-velocity 911 CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8) ) ! now k-velocity 912 CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8) ) ! now k-velocity 913 CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par ) ! now k-velocity 914 CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par ) ! now k-velocity 915 CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par ) ! now k-velocity 916 CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8) ) ! now k-velocity 917 CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8) ) ! now k-velocity 906 918 IF( ALLOCATED(ahtu) ) THEN 907 919 CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/dom_oce.F90
r10068 r11395 32 32 LOGICAL , PUBLIC :: ln_linssh !: =T linear free surface ==>> model level are fixed in time 33 33 LOGICAL , PUBLIC :: ln_meshmask !: =T create a mesh-mask file (mesh_mask.nc) 34 REAL(wp), PUBLIC :: rn_isfhmin !: threshold to discriminate grounded ice to floating ice35 34 REAL(wp), PUBLIC :: rn_rdt !: time step for the dynamics and tracer 36 35 REAL(wp), PUBLIC :: rn_atfp !: asselin time filter parameter … … 170 169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_h !: internal domain T-point mask (Figure 8.5 NEMO book) 171 170 172 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfdep !: top first ocean level (ISF) 173 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: top first wet T-, U-, V-, F-level (ISF) 174 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfdep !: Iceshelf draft (ISF) 171 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: top first wet T-, U-, V-, F-level (ISF) 175 172 176 173 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask !: surface mask at T-,U-, V- and F-pts … … 281 278 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(9) ) 282 279 ! 283 ALLOCATE( misfdep(jpi,jpj) , mikt(jpi,jpj) , miku(jpi,jpj) , & 284 & risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(10) ) 280 ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(10) ) 285 281 ! 286 282 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) , & -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domain.F90
r10425 r11395 140 140 ! Read in masks to define closed seas and lakes 141 141 ! 142 DO jj = 1, jpj ! depth of the iceshelves143 DO ji = 1, jpi144 ik = mikt(ji,jj)145 risfdep(ji,jj) = gdepw_0(ji,jj,ik)146 END DO147 END DO148 !149 142 ht_0(:,:) = 0._wp ! Reference ocean thickness 150 143 hu_0(:,:) = 0._wp … … 293 286 & nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, nn_euler , & 294 287 & ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios 295 NAMELIST/namdom/ ln_linssh, rn_ isfhmin, rn_rdt, rn_atfp, ln_crs, ln_meshmask288 NAMELIST/namdom/ ln_linssh, rn_rdt, rn_atfp, ln_crs, ln_meshmask 296 289 #if defined key_netcdf4 297 290 NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip … … 412 405 WRITE(numout,*) ' linear free surface (=T) ln_linssh = ', ln_linssh 413 406 WRITE(numout,*) ' create mesh/mask file ln_meshmask = ', ln_meshmask 414 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' [m]'415 407 WRITE(numout,*) ' ocean time step rn_rdt = ', rn_rdt 416 408 WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domwri.F90
r10425 r11395 16 16 !! dom_stiff : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate) 17 17 !!---------------------------------------------------------------------- 18 USE isf ! ice shelf 18 19 USE dom_oce ! ocean space and time domain 19 20 USE phycst , ONLY : rsmall … … 155 156 156 157 ! note that mbkt is set to 1 over land ==> use surface tmask 157 zprt(:,:) = ssmask(:,:) *REAL( mbkt(:,:) , wp )158 zprt(:,:) = REAL( mbkt(:,:) , wp ) 158 159 CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points 159 zprt(:,:) = ssmask(:,:) *REAL( mikt(:,:) , wp )160 zprt(:,:) = REAL( mikt(:,:) , wp ) 160 161 CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points 161 zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )162 CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 ) ! ! nb of ocean T-points163 162 ! ! vertical mesh 164 163 CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8 ) ! ! scale factors -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/divhor.F90
r10425 r11395 22 22 USE sbc_oce, ONLY : ln_rnf, ln_isf ! surface boundary condition: ocean 23 23 USE sbcrnf ! river runoff 24 USE sbcisf! ice shelf24 USE isfhdiv ! ice shelf 25 25 USE iscplhsb ! ice sheet / ocean coupling 26 26 USE iscplini ! ice sheet / ocean coupling … … 100 100 ! 101 101 #endif 102 IF( ln_isf ) CALL sbc_isf_div( hdivn ) !== ice shelf ==! (update hdivn field) 103 ! 104 IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdivn ) !== ice sheet ==! (update hdivn field) 102 IF( ln_isf ) CALL isf_hdiv( hdivn ) !== ice shelf ==! (update hdivn field) 105 103 ! 106 104 CALL lbc_lnk( 'divhor', hdivn, 'T', 1. ) ! (no sign change) -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynhpg.F90
r10491 r11395 31 31 !!---------------------------------------------------------------------- 32 32 USE oce ! ocean dynamics and tracers 33 USE isf ! ice shelf 33 34 USE sbc_oce ! surface variable (only for the flag with ice shelf) 34 35 USE dom_oce ! ocean space and time domain -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90
r10425 r11395 29 29 USE sbc_oce ! Surface boundary condition: ocean fields 30 30 USE sbcrnf ! river runoffs 31 USE sbcisf ! ice shelf31 USE isfnxt 32 32 USE phycst ! physical constants 33 33 USE dynadv ! dynamics: vector invariant versus flux form … … 241 241 ENDIF 242 242 END IF 243 244 IF ( ln_isf ) THEN ! if ice shelf melting 245 DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too 246 DO jj = 1, jpj 247 DO ji = 1, jpi 248 IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj) ) THEN 249 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 250 & * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk) 251 ELSEIF ( jk==misfkb(ji,jj) ) THEN 252 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 253 & * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk) 254 ENDIF 255 END DO 256 END DO 257 END DO 258 END IF 243 ! 244 ! ice shelf melting 245 IF ( ln_isf ) CALL isf_dynnxt( zcoef ) 259 246 ! 260 247 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90
r10742 r11395 33 33 USE zdf_oce ! vertical physics: variables 34 34 USE zdfdrg ! vertical physics: top/bottom drag coef. 35 USE sbcisf! ice shelf variable (fwfisf)35 USE isf ! ice shelf variable (fwfisf) 36 36 USE sbcapr ! surface boundary condition: atmospheric pressure 37 37 USE dynadv , ONLY: ln_dynadv_vec … … 632 632 ! ! Surface net water flux and rivers 633 633 IF (ln_bt_fw) THEN 634 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf (:,:) )634 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) 635 635 ELSE 636 636 zztmp = r1_rau0 * r1_2 637 637 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 638 & + fwfisf(:,:) + fwfisf_b(:,:) ) 638 & + fwfisf_cav(:,:) + fwfisf_cav_b(:,:) & 639 & + fwfisf_par(:,:) + fwfisf_par_b(:,:) ) 639 640 ENDIF 640 641 ! -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/sshwzv.F90
r11293 r11395 17 17 !!---------------------------------------------------------------------- 18 18 USE oce ! ocean dynamics and tracers variables 19 USE isf ! ice shelf 19 20 USE dom_oce ! ocean space and time domain variables 20 21 USE sbc_oce ! surface boundary condition: ocean … … 256 257 sshb(:,:) = sshb(:,:) - zcoef * ( emp_b(:,:) - emp (:,:) & 257 258 & - rnf_b(:,:) + rnf (:,:) & 258 & + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:) 259 & + fwfisf_cav_b(:,:) - fwfisf_cav(:,:) & 260 & + fwfisf_par_b(:,:) - fwfisf_par(:,:) ) * ssmask(:,:) 259 261 ENDIF 260 262 sshn(:,:) = ssha(:,:) ! now <-- after -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/iscplini.F90
r11310 r11395 23 23 PUBLIC iscpl_init 24 24 PUBLIC iscpl_alloc 25 26 25 ! !!* namsbc_iscpl namelist * 27 LOGICAL , PUBLIC :: ln_ hsb!:26 LOGICAL , PUBLIC :: ln_iscpl_hsb !: 28 27 INTEGER , PUBLIC :: nn_fiscpl !: 29 28 INTEGER , PUBLIC :: nn_drown !: 30 29 ! 31 30 INTEGER , PUBLIC :: nstp_iscpl !: 32 31 REAL(wp), PUBLIC :: rdt_iscpl !: … … 57 56 !!---------------------------------------------------------------------- 58 57 INTEGER :: ios ! Local integer output status for namelist read 59 NAMELIST/namsbc_iscpl/ nn_fiscpl, ln_ hsb, nn_drown58 NAMELIST/namsbc_iscpl/ nn_fiscpl, ln_iscpl_hsb, nn_drown 60 59 !!---------------------------------------------------------------------- 61 60 ! 62 61 nn_fiscpl = 0 63 ln_ hsb = .FALSE.62 ln_iscpl_hsb = .FALSE. 64 63 REWIND( numnam_ref ) ! Namelist namsbc_iscpl in reference namelist : Ice sheet coupling 65 64 READ ( numnam_ref, namsbc_iscpl, IOSTAT = ios, ERR = 901) … … 76 75 WRITE(numout,*) 'iscpl_rst:' 77 76 WRITE(numout,*) '~~~~~~~~~' 78 WRITE(numout,*) ' coupling flag (ln_iscpl ) = ', ln_iscpl79 WRITE(numout,*) ' conservation flag (ln_ hsb ) = ', ln_hsb80 WRITE(numout,*) ' nb of stp for cons (rn_fiscpl )= ', nstp_iscpl77 WRITE(numout,*) ' coupling flag (ln_iscpl ) = ', ln_iscpl 78 WRITE(numout,*) ' conservation flag (ln_iscpl_hsb) = ', ln_iscpl_hsb 79 WRITE(numout,*) ' nb of stp for cons (rn_fiscpl ) = ', nstp_iscpl 81 80 IF (nstp_iscpl .NE. nn_fiscpl) WRITE(numout,*) 'W A R N I N G: nb of stp for cons has been modified & 82 81 & (larger than run length)' -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/iscplrst.F90
r11310 r11395 71 71 CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b ) 72 72 ! 73 IF ( ln_ hsb ) THEN ! compute correction if conservation needed73 IF ( ln_iscpl_hsb ) THEN ! compute correction if conservation needed 74 74 IF( iscpl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' ) 75 75 CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl) … … 79 79 IF( ln_meshmask .AND. ln_iscpl ) CALL dom_wri 80 80 ! 81 IF ( ln_ hsb ) THEN81 IF ( ln_iscpl_hsb ) THEN 82 82 cfile='correction' 83 83 cfile = TRIM( cfile ) -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfmlt.F90
r11310 r11395 1 MODULE sbcisf1 MODULE isfmlt 2 2 !!====================================================================== 3 3 !! *** MODULE sbcisf *** 4 !! Surface module : update surface ocean boundary condition under ice 5 !! shelf 4 !! Surface module : compute iceshelf melt and heat flux 6 5 !!====================================================================== 7 6 !! History : 3.2 ! 2011-02 (C.Harris ) Original code isf cav 8 7 !! X.X ! 2006-02 (C. Wang ) Original code bg03 9 8 !! 3.4 ! 2013-03 (P. Mathiot) Merging + parametrization 10 !!---------------------------------------------------------------------- 11 12 !!---------------------------------------------------------------------- 13 !! sbc_isf : update sbc under ice shelf 9 !! 4.1 ! 2019-09 (P. Mathiot) Split param/explicit ice shelf and re-organisation 10 !!---------------------------------------------------------------------- 11 12 !!---------------------------------------------------------------------- 13 !! isfmlt : compute iceshelf melt and heat flux 14 14 !!---------------------------------------------------------------------- 15 15 USE oce ! ocean dynamics and tracers … … 25 25 USE lbclnk ! 26 26 USE lib_fortran ! glob_sum 27 ! 28 USE isfrst ! iceshelf restart 29 USE isftbl ! ice shelf boundary layer 30 USE isfpar ! ice shelf parametrisation 31 USE isfcav ! ice shelf cavity 32 USE isf ! isf variables 27 33 28 34 IMPLICIT NONE 35 29 36 PRIVATE 30 37 31 PUBLIC sbc_isf, sbc_isf_init, sbc_isf_div, sbc_isf_alloc ! routine called in sbcmod and divhor 32 33 ! public in order to be able to output then 34 35 REAL(wp), PUBLIC :: rn_hisf_tbl !: thickness of top boundary layer [m] 36 INTEGER , PUBLIC :: nn_isf !: flag to choose between explicit/param/specified 37 INTEGER , PUBLIC :: nn_isfblk !: flag to choose the bulk formulation to compute the ice shelf melting 38 INTEGER , PUBLIC :: nn_gammablk !: flag to choose how the exchange coefficient is computed 39 REAL(wp), PUBLIC :: rn_gammat0 !: temperature exchange coeficient [] 40 REAL(wp), PUBLIC :: rn_gammas0 !: salinity exchange coeficient [] 41 42 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt , misfkb !: Level of ice shelf base 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rzisf_tbl !: depth of calving front (shallowest point) nn_isf ==2/3 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl, rhisf_tbl_0 !: thickness of tbl [m] 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hisf_tbl !: 1/thickness of tbl 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ralpha !: proportion of bottom cell influenced by tbl 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfLeff !: effective length (Leff) BG03 nn_isf==2 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ttbl, stbl, utbl, vtbl !: top boundary layer variable at T point 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf [W/m2] 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_tsc_b, risf_tsc !: before and now T & S isf contents [K.m/s & PSU.m/s] 51 52 LOGICAL, PUBLIC :: l_isfcpl = .false. !: isf recieved from oasis 53 54 REAL(wp), PUBLIC, SAVE :: rcpisf = 2000.0_wp !: specific heat of ice shelf [J/kg/K] 55 REAL(wp), PUBLIC, SAVE :: rkappa = 1.54e-6_wp !: heat diffusivity through the ice-shelf [m2/s] 56 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp !: volumic mass of ice shelf [kg/m3] 57 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp !: air temperature on top of ice shelf [C] 58 REAL(wp), PUBLIC, SAVE :: rLfusisf = 0.334e6_wp !: latent heat of fusion of ice shelf [J/kg] 59 60 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 61 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 62 TYPE(FLD_N) , PUBLIC :: sn_fwfisf !: information about the isf melting file to be read 63 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 64 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the isf melting param. file to be read 65 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf 66 TYPE(FLD_N) , PUBLIC :: sn_depmax_isf !: information about the grounding line depth file to be read 67 TYPE(FLD_N) , PUBLIC :: sn_depmin_isf !: information about the calving line depth file to be read 68 TYPE(FLD_N) , PUBLIC :: sn_Leff_isf !: information about the effective length file to be read 69 38 PUBLIC isf_mlt, isf_mlt_init ! routine called in sbcmod and divhor 39 70 40 !!---------------------------------------------------------------------- 71 41 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 75 45 CONTAINS 76 46 77 SUBROUTINE sbc_isf( kt )47 SUBROUTINE isf_mlt( kt ) 78 48 !!--------------------------------------------------------------------- 79 !! *** ROUTINE sbc_isf *** 80 !! 81 !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 82 !! melting and freezing 83 !! 84 !! ** Method : 4 parameterizations are available according to nn_isf 85 !! nn_isf = 1 : Realistic ice_shelf formulation 86 !! 2 : Beckmann & Goose parameterization 87 !! 3 : Specified runoff in deptht (Mathiot & al. ) 88 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 49 !! *** ROUTINE isf_mlt *** 50 !! 51 !! ** Purpose : 52 !! 53 !! ** Method : 54 !! 89 55 !!---------------------------------------------------------------------- 90 56 INTEGER, INTENT(in) :: kt ! ocean time step 91 !92 INTEGER :: ji, jj, jk ! loop index93 INTEGER :: ikt, ikb ! local integers94 REAL(wp), DIMENSION(jpi,jpj) :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)95 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zqhcisf2d96 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zfwfisf3d, zqhcisf3d, zqlatisf3d97 57 !!--------------------------------------------------------------------- 98 58 ! 99 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN ! compute salt and heat flux 100 ! 101 SELECT CASE ( nn_isf ) 102 CASE ( 1 ) ! realistic ice shelf formulation 103 ! compute T/S/U/V for the top boundary layer 104 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 105 CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 106 CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U') 107 CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V') 108 ! iom print 109 CALL iom_put('ttbl',ttbl(:,:)) 110 CALL iom_put('stbl',stbl(:,:)) 111 CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 112 CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 113 ! compute fwf and heat flux 114 ! compute fwf and heat flux 115 IF( .NOT.l_isfcpl ) THEN ; CALL sbc_isf_cav (kt) 116 ELSE ; qisf(:,:) = fwfisf(:,:) * rLfusisf ! heat flux 117 ENDIF 118 ! 119 CASE ( 2 ) ! Beckmann and Goosse parametrisation 120 stbl(:,:) = soce 121 CALL sbc_isf_bg03(kt) 122 ! 123 CASE ( 3 ) ! specified runoff in depth (Mathiot et al., XXXX in preparation) 124 ! specified runoff in depth (Mathiot et al., XXXX in preparation) 125 IF( .NOT.l_isfcpl ) THEN 126 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 127 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) 128 ENDIF 129 qisf(:,:) = fwfisf(:,:) * rLfusisf ! heat flux 130 stbl(:,:) = soce 131 ! 132 CASE ( 4 ) ! specified fwf and heat flux forcing beneath the ice shelf 133 ! ! specified fwf and heat flux forcing beneath the ice shelf 134 IF( .NOT.l_isfcpl ) THEN 135 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 136 !CALL fld_read ( kt, nn_fsbc, sf_qisf ) 137 fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1) ! fwf 138 ENDIF 139 qisf(:,:) = fwfisf(:,:) * rLfusisf ! heat flux 140 stbl(:,:) = soce 141 ! 142 END SELECT 143 144 ! compute tsc due to isf 145 ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU. 146 ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0). 147 ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3) 148 DO jj = 1,jpj 149 DO ji = 1,jpi 150 zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj)) 151 END DO 152 END DO 153 CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 154 155 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 156 risf_tsc(:,:,jp_sal) = 0.0_wp 157 158 ! lbclnk 159 CALL lbc_lnk_multi( 'sbcisf', risf_tsc(:,:,jp_tem), 'T', 1., risf_tsc(:,:,jp_sal), 'T', 1., fwfisf,'T', 1., qisf, 'T', 1.) 160 ! output 161 IF( iom_use('iceshelf_cea') ) CALL iom_put( 'iceshelf_cea', -fwfisf(:,:) ) ! isf mass flux 162 IF( iom_use('hflx_isf_cea') ) CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * rau0 * rcp ) ! isf sensible+latent heat (W/m2) 163 IF( iom_use('qlatisf' ) ) CALL iom_put( 'qlatisf' , qisf(:,:) ) ! isf latent heat 164 IF( iom_use('fwfisf' ) ) CALL iom_put( 'fwfisf' , fwfisf(:,:) ) ! isf mass flux (opposite sign) 165 166 ! Diagnostics 167 IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 168 ALLOCATE( zfwfisf3d(jpi,jpj,jpk) , zqhcisf3d(jpi,jpj,jpk) , zqlatisf3d(jpi,jpj,jpk) ) 169 ALLOCATE( zqhcisf2d(jpi,jpj) ) 170 ! 171 zfwfisf3d (:,:,:) = 0._wp ! 3d ice shelf melting (kg/m2/s) 172 zqhcisf3d (:,:,:) = 0._wp ! 3d heat content flux (W/m2) 173 zqlatisf3d(:,:,:) = 0._wp ! 3d ice shelf melting latent heat flux (W/m2) 174 zqhcisf2d (:,:) = fwfisf(:,:) * zt_frz * rcp ! 2d heat content flux (W/m2) 175 ! 176 DO jj = 1,jpj 177 DO ji = 1,jpi 178 ikt = misfkt(ji,jj) 179 ikb = misfkb(ji,jj) 180 DO jk = ikt, ikb - 1 181 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 182 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 183 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 184 END DO 185 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) & 186 & * ralpha(ji,jj) * e3t_n(ji,jj,jk) 187 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) & 188 & * ralpha(ji,jj) * e3t_n(ji,jj,jk) 189 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) & 190 & * ralpha(ji,jj) * e3t_n(ji,jj,jk) 191 END DO 192 END DO 193 ! 194 CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:)) 195 CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:)) 196 CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:)) 197 CALL iom_put('qhcisf' , zqhcisf2d (:,: )) 198 ! 199 DEALLOCATE( zfwfisf3d, zqhcisf3d, zqlatisf3d ) 200 DEALLOCATE( zqhcisf2d ) 201 ENDIF 202 ! 203 ENDIF 204 205 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! 206 IF( ln_rstart .AND. & ! Restart: read in restart file 207 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN 208 IF(lwp) WRITE(numout,*) ' nit000-1 isf tracer content forcing fields read in the restart file' 209 CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) , ldxios = lrxios ) ! before salt content isf_tsc trend 210 CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b' , risf_tsc_b(:,:,jp_sal), ldxios = lrxios ) ! before salt content isf_tsc trend 211 CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b' , risf_tsc_b(:,:,jp_tem), ldxios = lrxios ) ! before salt content isf_tsc trend 212 ELSE 213 fwfisf_b(:,:) = fwfisf(:,:) 214 risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 215 ENDIF 216 ENDIF 59 IF ( ln_isfcav_mlt ) THEN 60 ! 61 ! before time step 62 IF ( kt /= nit000 ) THEN 63 risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:) 64 fwfisf_cav_b(:,:) = fwfisf_cav(:,:) 65 END IF 66 ! 67 ! compute tbl lvl/h 68 CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav) 69 ! 70 ! compute ice shelf melt 71 CALL isf_cav( kt, risf_cav_tsc, fwfisf_cav) 72 ! 73 ! write restart variables (qoceisf, qhcisf, fwfisf for now and before) 74 IF (lrst_oce) CALL isfrst_write(kt, 'cav', risf_cav_tsc, fwfisf_cav) 75 ! 76 END IF 217 77 ! 218 IF( lrst_oce ) THEN 219 IF(lwp) WRITE(numout,*) 220 IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ', & 221 & 'at it= ', kt,' date= ', ndastp 222 IF(lwp) WRITE(numout,*) '~~~~' 223 IF( lwxios ) CALL iom_swap( cwxios_context ) 224 CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) , ldxios = lwxios ) 225 CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem), ldxios = lwxios ) 226 CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal), ldxios = lwxios ) 227 IF( lwxios ) CALL iom_swap( cxios_context ) 228 ENDIF 229 ! 230 END SUBROUTINE sbc_isf 231 232 233 INTEGER FUNCTION sbc_isf_alloc() 234 !!---------------------------------------------------------------------- 235 !! *** FUNCTION sbc_isf_rnf_alloc *** 236 !!---------------------------------------------------------------------- 237 sbc_isf_alloc = 0 ! set to zero if no array to be allocated 238 IF( .NOT. ALLOCATED( qisf ) ) THEN 239 ALLOCATE( risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj) , & 240 & rhisf_tbl(jpi,jpj) , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj) , & 241 & ttbl(jpi,jpj) , stbl(jpi,jpj) , utbl(jpi,jpj) , & 242 & vtbl(jpi, jpj) , risfLeff(jpi,jpj) , rhisf_tbl_0(jpi,jpj), & 243 & ralpha(jpi,jpj) , misfkt(jpi,jpj) , misfkb(jpi,jpj) , & 244 & STAT= sbc_isf_alloc ) 245 ! 246 CALL mpp_sum ( 'sbcisf', sbc_isf_alloc ) 247 IF( sbc_isf_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf_alloc: failed to allocate arrays.' ) 248 ! 249 ENDIF 250 END FUNCTION 251 252 253 SUBROUTINE sbc_isf_init 78 IF ( ln_isfpar_mlt ) THEN 79 ! 80 ! before time step 81 IF ( kt /= nit000 ) THEN 82 risf_par_tsc_b (:,:) = risf_par_tsc (:,:) 83 fwfisf_par_b(:,:) = fwfisf_par(:,:) 84 END IF 85 ! 86 ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl) 87 CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par) 88 ! 89 ! compute ice shelf melt 90 CALL isf_par( kt, risf_par_tsc, fwfisf_par) 91 ! 92 ! write restart variables (qoceisf, qhcisf, fwfisf for now and before) 93 IF (lrst_oce) CALL isfrst_write(kt, 'par', risf_par_tsc, fwfisf_par) 94 ! 95 END IF 96 ! 97 END SUBROUTINE isf_mlt 98 99 SUBROUTINE isf_mlt_init 254 100 !!--------------------------------------------------------------------- 255 !! *** ROUTINE sbc_isf_init *** 256 !! 257 !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation 258 !! 259 !! ** Method : 4 parameterizations are available according to nn_isf 260 !! nn_isf = 1 : Realistic ice_shelf formulation 261 !! 2 : Beckmann & Goose parameterization 262 !! 3 : Specified runoff in deptht (Mathiot & al. ) 263 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 264 !!---------------------------------------------------------------------- 265 INTEGER :: ji, jj, jk ! loop index 266 INTEGER :: ik ! current level index 267 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 101 !! *** ROUTINE isfmlt_init *** 102 !! 103 !! ** Purpose : 104 !! 105 !! ** Method : 106 !! 107 !!---------------------------------------------------------------------- 268 108 INTEGER :: inum, ierror 269 109 INTEGER :: ios ! Local integer output status for namelist read 270 REAL(wp) :: zhk 271 CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file 272 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale 273 !!---------------------------------------------------------------------- 274 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, & 275 & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 110 INTEGER :: ikt, ikb 111 INTEGER :: ji, jj 112 !!---------------------------------------------------------------------- 113 NAMELIST/namisf/ ln_isfcav_mlt, cn_isfcav_mlt, cn_gammablk, rn_gammat0, rn_gammas0, rn_htbl, sn_isfcav_fwf, & 114 & ln_isfpar_mlt, cn_isfpar_mlt, sn_isfpar_fwf, sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff 276 115 !!---------------------------------------------------------------------- 277 116 278 117 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 279 READ ( numnam_ref, nam sbc_isf, IOSTAT = ios, ERR = 901)280 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam sbc_isf in reference namelist', lwp )118 READ ( numnam_ref, namisf, IOSTAT = ios, ERR = 901) 119 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namisf in reference namelist', lwp ) 281 120 282 121 REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs 283 READ ( numnam_cfg, nam sbc_isf, IOSTAT = ios, ERR = 902 )284 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam sbc_isf in configuration namelist', lwp )285 IF(lwm) WRITE ( numond, nam sbc_isf )286 122 READ ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 ) 123 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namisf in configuration namelist', lwp ) 124 IF(lwm) WRITE ( numond, namisf ) 125 ! 287 126 IF(lwp) WRITE(numout,*) 288 IF(lwp) WRITE(numout,*) ' sbc_isf_init : heat flux of the ice shelf'127 IF(lwp) WRITE(numout,*) 'isf_init : ice shelf initialisation' 289 128 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 290 IF(lwp) WRITE(numout,*) ' Namelist namsbc_isf :' 291 IF(lwp) WRITE(numout,*) ' type ice shelf melting/freezing nn_isf = ', nn_isf 292 IF(lwp) WRITE(numout,*) ' bulk formulation (nn_isf=1 only) nn_isfblk = ', nn_isfblk 293 IF(lwp) WRITE(numout,*) ' thickness of the top boundary layer rn_hisf_tbl = ', rn_hisf_tbl 294 IF(lwp) WRITE(numout,*) ' gamma formulation nn_gammablk = ', nn_gammablk 295 IF(lwp) WRITE(numout,*) ' gammat coefficient rn_gammat0 = ', rn_gammat0 296 IF(lwp) WRITE(numout,*) ' gammas coefficient rn_gammas0 = ', rn_gammas0 297 IF(lwp) WRITE(numout,*) ' top drag coef. used (from namdrg_top) rn_Cd0 = ', r_Cdmin_top 298 299 300 ! 1 = presence of ISF 2 = bg03 parametrisation 301 ! 3 = rnf file for isf 4 = ISF fwf specified 302 ! option 1 and 4 need ln_isfcav = .true. (domzgr) 303 ! 304 ! Allocate public variable 305 IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 129 IF(lwp) WRITE(numout,*) ' Namelist namisf :' 130 IF(lwp) WRITE(numout,*) ' melt inside the cavity ln_isfcav_mlt = ', ln_isfcav_mlt 131 IF ( ln_isfcav ) THEN 132 IF(lwp) WRITE(numout,*) ' melt formulation cn_isfcav_mlt = ', cn_isfcav_mlt 133 IF(lwp) WRITE(numout,*) ' thickness of the top boundary layer rn_htbl = ', rn_htbl 134 IF(lwp) WRITE(numout,*) ' gamma formulation cn_gammablk = ', cn_gammablk 135 IF ( TRIM(cn_gammablk) .NE. 'spe' ) THEN 136 IF(lwp) WRITE(numout,*) ' gammat coefficient rn_gammat0 = ', rn_gammat0 137 IF(lwp) WRITE(numout,*) ' gammas coefficient rn_gammas0 = ', rn_gammas0 138 IF(lwp) WRITE(numout,*) ' top drag coef. used (from namdrg_top) rn_Cd0 = ', r_Cdmin_top 139 END IF 140 END IF 141 IF(lwp) WRITE(numout,*) '' 142 IF(lwp) WRITE(numout,*) ' ice shelf melt parametrisation ln_isfpar_mlt = ', ln_isfpar_mlt 143 IF ( ln_isfpar_mlt ) THEN 144 IF(lwp) WRITE(numout,*) ' isf parametrisation formulation cn_isfpar_mlt = ', cn_isfpar_mlt 145 END IF 146 IF(lwp) WRITE(numout,*) '' 147 ! 148 ! Allocate public array 149 CALL isf_alloc() 306 150 ! 307 151 ! initialisation 308 qisf (:,:) = 0._wp ; fwfisf (:,:) = 0._wp 309 risf_tsc(:,:,:) = 0._wp ; fwfisf_b(:,:) = 0._wp 310 ! 311 ! define isf tbl tickness, top and bottom indice 312 SELECT CASE ( nn_isf ) 313 CASE ( 1 ) 314 IF(lwp) WRITE(numout,*) 315 IF(lwp) WRITE(numout,*) ' ==>>> presence of under iceshelf seas (nn_isf = 1)' 316 rhisf_tbl(:,:) = rn_hisf_tbl 317 misfkt (:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 318 ! 319 CASE ( 2 , 3 ) 320 IF( .NOT.l_isfcpl ) THEN 321 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 322 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 323 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 324 ENDIF 325 ! read effective lenght (BG03) 326 IF( nn_isf == 2 ) THEN 327 IF(lwp) WRITE(numout,*) 328 IF(lwp) WRITE(numout,*) ' ==>>> bg03 parametrisation (nn_isf = 2)' 329 CALL iom_open( sn_Leff_isf%clname, inum ) 330 cvarLeff = TRIM(sn_Leff_isf%clvar) 331 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 332 CALL iom_close(inum) 333 ! 334 risfLeff = risfLeff*1000.0_wp !: convertion in m 335 ELSE 336 IF(lwp) WRITE(numout,*) 337 IF(lwp) WRITE(numout,*) ' ==>>> rnf file for isf (nn_isf = 3)' 338 ENDIF 339 ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 340 CALL iom_open( sn_depmax_isf%clname, inum ) 341 cvarhisf = TRIM(sn_depmax_isf%clvar) 342 CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 343 CALL iom_close(inum) 344 ! 345 CALL iom_open( sn_depmin_isf%clname, inum ) 346 cvarzisf = TRIM(sn_depmin_isf%clvar) 347 CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 348 CALL iom_close(inum) 349 ! 350 rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer 351 352 !! compute first level of the top boundary layer 353 DO ji = 1, jpi 354 DO jj = 1, jpj 355 ik = 2 356 !!gm potential bug: use gdepw_0 not _n 357 DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ; ik = ik + 1 ; END DO 358 misfkt(ji,jj) = ik-1 359 END DO 360 END DO 361 ! 362 CASE ( 4 ) 363 IF(lwp) WRITE(numout,*) 364 IF(lwp) WRITE(numout,*) ' ==>>> specified fresh water flux in ISF (nn_isf = 4)' 365 ! as in nn_isf == 1 366 rhisf_tbl(:,:) = rn_hisf_tbl 367 misfkt (:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 368 ! 369 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 370 IF( .NOT.l_isfcpl ) THEN 371 ALLOCATE( sf_fwfisf(1), STAT=ierror ) 372 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 373 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 374 ENDIF 375 ! 376 CASE DEFAULT 377 CALL ctl_stop( 'sbc_isf_init: wrong value of nn_isf' ) 378 END SELECT 379 380 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 381 382 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 152 ! cav 153 misfkt_cav(:,:) = mikt(:,:) ; misfkb_cav(:,:) = mikt(:,:) 154 fwfisf_cav(:,:) = 0.0_wp ; fwfisf_cav_b(:,:) = 0.0_wp 155 rhisf_tbl_cav(:,:) = 1e-20 ; rfrac_tbl_cav(:,:) = 0.0_wp 156 risf_cav_tsc(:,:,:) = 0.0_wp ; risf_cav_tsc_b(:,:,:) = 0.0_wp 157 ! 158 mskisf_cav(:,:) = (1._wp - tmask(:,:,1)) * ssmask(:,:) 159 ! 160 ! par 161 misfkt_par(:,:) = 1 ; misfkb_par(:,:) = 1 162 fwfisf_par(:,:) = 0.0_wp ; fwfisf_par_b(:,:) = 0.0_wp 163 rhisf_tbl_par(:,:) = 1e-20 ; rfrac_tbl_par(:,:) = 0.0_wp 164 risf_par_tsc(:,:,:) = 0.0_wp ; risf_par_tsc_b(:,:,:) = 0.0_wp 165 ! 166 mskisf_par(:,:) = 0 167 ! 168 ! initialisation useful variable 169 r1_Lfusisf = 1._wp / rLfusisf 170 ! 383 171 DO jj = 1,jpj 384 172 DO ji = 1,jpi 385 ikt = misfkt(ji,jj) 386 ikb = misfkt(ji,jj) 387 ! thickness of boundary layer at least the top level thickness 388 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt)) 389 390 ! determine the deepest level influenced by the boundary layer 391 DO jk = ikt+1, mbkt(ji,jj) 392 IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 393 END DO 394 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 395 misfkb(ji,jj) = ikb ! last wet level of the tbl 396 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 397 398 zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 399 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 173 ikt = mikt(ji,jj) 174 ikb = mbkt(ji,jj) 175 bathy (ji,jj) = gdepw_0(ji,jj,ikb+1) 176 risfdep(ji,jj) = gdepw_0(ji,jj,ikt ) 400 177 END DO 401 178 END DO 402 403 IF( lwxios ) THEN 404 CALL iom_set_rstw_var_active('fwf_isf_b') 405 CALL iom_set_rstw_var_active('isf_hc_b') 406 CALL iom_set_rstw_var_active('isf_sc_b') 407 ENDIF 408 409 410 END SUBROUTINE sbc_isf_init 411 412 413 SUBROUTINE sbc_isf_bg03(kt) 414 !!--------------------------------------------------------------------- 415 !! *** ROUTINE sbc_isf_bg03 *** 416 !! 417 !! ** Purpose : add net heat and fresh water flux from ice shelf melting 418 !! into the adjacent ocean 419 !! 420 !! ** Method : See reference 421 !! 422 !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 423 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 424 !! (hereafter BG) 425 !! History : 06-02 (C. Wang) Original code 426 !!---------------------------------------------------------------------- 427 INTEGER, INTENT ( in ) :: kt 428 ! 429 INTEGER :: ji, jj, jk ! dummy loop index 430 INTEGER :: ik ! current level 431 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 432 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 433 REAL(wp) :: zt_frz ! freezing point temperature at depth z 434 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 435 !!---------------------------------------------------------------------- 436 ! 437 DO ji = 1, jpi 438 DO jj = 1, jpj 439 ik = misfkt(ji,jj) 440 !! Initialize arrays to 0 (each step) 441 zt_sum = 0.e0_wp 442 IF ( ik > 1 ) THEN 443 ! 1. -----------the average temperature between 200m and 600m --------------------- 444 DO jk = misfkt(ji,jj),misfkb(ji,jj) 445 ! Calculate freezing temperature 446 zpress = grav*rau0*gdept_n(ji,jj,ik)*1.e-04 447 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 448 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp 449 END DO 450 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 451 ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 452 ! For those corresponding to zonal boundary 453 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 454 & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) 455 456 fwfisf(ji,jj) = qisf(ji,jj) / rLfusisf !fresh water flux kg/(m2s) 457 fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 458 !add to salinity trend 459 ELSE 460 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 461 END IF 462 END DO 463 END DO 464 ! 465 END SUBROUTINE sbc_isf_bg03 466 467 468 SUBROUTINE sbc_isf_cav( kt ) 469 !!--------------------------------------------------------------------- 470 !! *** ROUTINE sbc_isf_cav *** 471 !! 472 !! ** Purpose : handle surface boundary condition under ice shelf 473 !! 474 !! ** Method : - 475 !! 476 !! ** Action : utau, vtau : remain unchanged 477 !! taum, wndm : remain unchanged 478 !! qns : update heat flux below ice shelf 479 !! emp, emps : update freshwater flux below ice shelf 480 !!--------------------------------------------------------------------- 481 INTEGER, INTENT(in) :: kt ! ocean time step 482 ! 483 INTEGER :: ji, jj ! dummy loop indices 484 INTEGER :: nit 485 LOGICAL :: lit 486 REAL(wp) :: zlamb1, zlamb2, zlamb3 487 REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 488 REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac 489 REAL(wp) :: zeps = 1.e-20_wp 490 REAL(wp) :: zerr 491 REAL(wp), DIMENSION(jpi,jpj) :: zfrz 492 REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas 493 REAL(wp), DIMENSION(jpi,jpj) :: zfwflx, zhtflx, zhtflx_b 494 !!--------------------------------------------------------------------- 495 ! 496 ! coeficient for linearisation of potential tfreez 497 ! Crude approximation for pressure (but commonly used) 498 IF ( l_useCT ) THEN ! linearisation from Jourdain et al. (2017) 499 zlamb1 =-0.0564_wp 500 zlamb2 = 0.0773_wp 501 zlamb3 =-7.8633e-8 * grav * rau0 502 ELSE ! linearisation from table 4 (Asay-Davis et al., 2015) 503 zlamb1 =-0.0573_wp 504 zlamb2 = 0.0832_wp 505 zlamb3 =-7.53e-8 * grav * rau0 506 ENDIF 507 ! 508 ! initialisation 509 zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 510 zhtflx (:,:) = 0.0_wp ; zhtflx_b(:,:) = 0.0_wp 511 zfwflx (:,:) = 0.0_wp 512 513 ! compute ice shelf melting 514 nit = 1 ; lit = .TRUE. 515 DO WHILE ( lit ) ! maybe just a constant number of iteration as in blk_core is fine 516 SELECT CASE ( nn_isfblk ) 517 CASE ( 1 ) ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 518 ! Calculate freezing temperature 519 CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) 520 521 ! compute gammat every where (2d) 522 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 523 524 ! compute upward heat flux zhtflx and upward water flux zwflx 525 DO jj = 1, jpj 526 DO ji = 1, jpi 527 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 528 zfwflx(ji,jj) = - zhtflx(ji,jj)/rLfusisf 529 END DO 530 END DO 531 532 ! Compute heat flux and upward fresh water flux 533 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 534 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 535 536 CASE ( 2 ) ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 537 ! compute gammat every where (2d) 538 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 539 540 ! compute upward heat flux zhtflx and upward water flux zwflx 541 ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015) 542 DO jj = 1, jpj 543 DO ji = 1, jpi 544 ! compute coeficient to solve the 2nd order equation 545 zeps1 = rcp*rau0*zgammat(ji,jj) 546 zeps2 = rLfusisf*rau0*zgammas(ji,jj) 547 zeps3 = rhoisf*rcpisf*rkappa/MAX(risfdep(ji,jj),zeps) 548 zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 549 zeps6 = zeps4-ttbl(ji,jj) 550 zeps7 = zeps4-tsurf 551 zaqe = zlamb1 * (zeps1 + zeps3) 552 zaqer = 0.5_wp/MIN(zaqe,-zeps) 553 zbqe = zeps1*zeps6+zeps3*zeps7-zeps2 554 zcqe = zeps2*stbl(ji,jj) 555 zdis = zbqe*zbqe-4.0_wp*zaqe*zcqe 556 557 ! Presumably zdis can never be negative because gammas is very small compared to gammat 558 ! compute s freeze 559 zsfrz=(-zbqe-SQRT(zdis))*zaqer 560 IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 561 562 ! compute t freeze (eq. 22) 563 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 564 565 ! zfwflx is upward water flux 566 ! zhtflx is upward heat flux (out of ocean) 567 ! compute the upward water and heat flux (eq. 28 and eq. 29) 568 zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) 569 zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 570 END DO 571 END DO 572 573 ! compute heat and water flux 574 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 575 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 576 577 END SELECT 578 579 ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 580 IF ( nn_gammablk < 2 ) THEN ; lit = .FALSE. 581 ELSE 582 ! check total number of iteration 583 IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 584 ELSE ; nit = nit + 1 585 END IF 586 587 ! compute error between 2 iterations 588 ! if needed save gammat and compute zhtflx_b for next iteration 589 zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 590 IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 591 ELSE ; zhtflx_b(:,:) = zhtflx(:,:) 592 END IF 593 END IF 594 END DO 595 ! 596 CALL iom_put('isfgammat', zgammat) 597 CALL iom_put('isfgammas', zgammas) 598 ! 599 END SUBROUTINE sbc_isf_cav 600 601 602 SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 603 !!---------------------------------------------------------------------- 604 !! ** Purpose : compute the coefficient echange for heat flux 605 !! 606 !! ** Method : gamma assume constant or depends of u* and stability 607 !! 608 !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 609 !! Jenkins et al., 2010, JPO, p2298-2312 610 !!--------------------------------------------------------------------- 611 REAL(wp), DIMENSION(:,:), INTENT( out) :: pgt , pgs ! 612 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf ! 613 ! 614 INTEGER :: ji, jj ! loop index 615 INTEGER :: ikt ! local integer 616 REAL(wp) :: zdku, zdkv ! U, V shear 617 REAL(wp) :: zPr, zSc, zRc ! Prandtl, Scmidth and Richardson number 618 REAL(wp) :: zmob, zmols ! Monin Obukov length, coriolis factor at T point 619 REAL(wp) :: zbuofdep, zhnu ! Bouyancy length scale, sublayer tickness 620 REAL(wp) :: zhmax ! limitation of mol 621 REAL(wp) :: zetastar ! stability parameter 622 REAL(wp) :: zgmolet, zgmoles, zgturb ! contribution of modelecular sublayer and turbulence 623 REAL(wp) :: zcoef ! temporary coef 624 REAL(wp) :: zdep 625 REAL(wp) :: zeps = 1.0e-20_wp 626 REAL(wp), PARAMETER :: zxsiN = 0.052_wp ! dimensionless constant 627 REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 628 REAL(wp), DIMENSION(2) :: zts, zab 629 REAL(wp), DIMENSION(jpi,jpj) :: zustar ! U, V at T point and friction velocity 630 !!--------------------------------------------------------------------- 631 ! 632 SELECT CASE ( nn_gammablk ) 633 CASE ( 0 ) ! gamma is constant (specified in namelist) 634 !! ISOMIP formulation (Hunter et al, 2006) 635 pgt(:,:) = rn_gammat0 636 pgs(:,:) = rn_gammas0 637 638 CASE ( 1 ) ! gamma is assume to be proportional to u* 639 !! Jenkins et al., 2010, JPO, p2298-2312 640 !! Adopted by Asay-Davis et al. (2015) 641 !! compute ustar (eq. 24) 642 !!gm NB use pCdU here so that it will incorporate local boost of Cd0 and log layer case : 643 !! zustar(:,:) = SQRT( rCdU_top(:,:) * SQRT(utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 644 !! or better : compute ustar in zdfdrg and use it here as well as in TKE, GLS and Co 645 !! 646 !! ===>>>> GM to be done this chrismas 647 !! 648 !!gm end 649 zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 650 651 !! Compute gammats 652 pgt(:,:) = zustar(:,:) * rn_gammat0 653 pgs(:,:) = zustar(:,:) * rn_gammas0 654 655 CASE ( 2 ) ! gamma depends of stability of boundary layer 656 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 657 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 658 !! compute ustar 659 zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 660 661 !! compute Pr and Sc number (can be improved) 662 zPr = 13.8_wp 663 zSc = 2432.0_wp 664 665 !! compute gamma mole 666 zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 667 zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 668 669 !! compute gamma 670 DO ji = 2, jpi 671 DO jj = 2, jpj 672 ikt = mikt(ji,jj) 673 674 IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think 675 pgt = rn_gammat0 676 pgs = rn_gammas0 677 ELSE 678 !! compute Rc number (as done in zdfric.F90) 679 !!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation 680 !!gm moreover, use Max(rn2,0) to take care of static instabilities.... 681 zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1) 682 ! ! shear of horizontal velocity 683 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & 684 & -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) ) 685 zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) & 686 & -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) ) 687 ! ! richardson number (minimum value set to zero) 688 zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 689 690 !! compute bouyancy 691 zts(jp_tem) = ttbl(ji,jj) 692 zts(jp_sal) = stbl(ji,jj) 693 zdep = gdepw_n(ji,jj,ikt) 694 ! 695 CALL eos_rab( zts, zdep, zab ) 696 ! 697 !! compute length scale 698 zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!! 699 700 !! compute Monin Obukov Length 701 ! Maximum boundary layer depth 702 zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp 703 ! Compute Monin obukhov length scale at the surface and Ekman depth: 704 zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 705 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 706 707 !! compute eta* (stability parameter) 708 zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp))) 709 710 !! compute the sublayer thickness 711 zhnu = 5 * znu / zustar(ji,jj) 712 713 !! compute gamma turb 714 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) & 715 & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 716 717 !! compute gammats 718 pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 719 pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 720 END IF 721 END DO 722 END DO 723 CALL lbc_lnk_multi( 'sbcisf', pgt, 'T', 1., pgs, 'T', 1.) 724 END SELECT 725 ! 726 END SUBROUTINE sbc_isf_gammats 727 728 729 SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 730 !!---------------------------------------------------------------------- 731 !! *** SUBROUTINE sbc_isf_tbl *** 732 !! 733 !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point 734 !! 735 !!---------------------------------------------------------------------- 736 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pvarin 737 REAL(wp), DIMENSION(:,:) , INTENT( out) :: pvarout 738 CHARACTER(len=1), INTENT(in ) :: cd_ptin ! point of variable in/out 739 ! 740 INTEGER :: ji, jj, jk ! loop index 741 INTEGER :: ikt, ikb ! top and bottom index of the tbl 742 REAL(wp) :: ze3, zhk 743 REAL(wp), DIMENSION(jpi,jpj) :: zhisf_tbl ! thickness of the tbl 744 !!---------------------------------------------------------------------- 745 746 ! initialisation 747 pvarout(:,:)=0._wp 748 749 SELECT CASE ( cd_ptin ) 750 CASE ( 'U' ) ! compute U in the top boundary layer at T- point 751 DO jj = 1,jpj 752 DO ji = 1,jpi 753 ikt = miku(ji,jj) ; ikb = miku(ji,jj) 754 ! thickness of boundary layer at least the top level thickness 755 zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) ) 756 757 ! determine the deepest level influenced by the boundary layer 758 DO jk = ikt+1, mbku(ji,jj) 759 IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 760 END DO 761 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 762 763 ! level fully include in the ice shelf boundary layer 764 DO jk = ikt, ikb - 1 765 ze3 = e3u_n(ji,jj,jk) 766 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 767 END DO 768 769 ! level partially include in ice shelf boundary layer 770 zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 771 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 772 END DO 773 END DO 774 DO jj = 2, jpj 775 DO ji = 2, jpi 776 !!gm a wet-point only average should be used here !!! 777 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 778 END DO 779 END DO 780 CALL lbc_lnk('sbcisf', pvarout,'T',-1.) 781 782 CASE ( 'V' ) ! compute V in the top boundary layer at T- point 783 DO jj = 1,jpj 784 DO ji = 1,jpi 785 ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 786 ! thickness of boundary layer at least the top level thickness 787 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt)) 788 789 ! determine the deepest level influenced by the boundary layer 790 DO jk = ikt+1, mbkv(ji,jj) 791 IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 792 END DO 793 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 794 795 ! level fully include in the ice shelf boundary layer 796 DO jk = ikt, ikb - 1 797 ze3 = e3v_n(ji,jj,jk) 798 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 799 END DO 800 801 ! level partially include in ice shelf boundary layer 802 zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 803 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 804 END DO 805 END DO 806 DO jj = 2, jpj 807 DO ji = 2, jpi 808 !!gm a wet-point only average should be used here !!! 809 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 810 END DO 811 END DO 812 CALL lbc_lnk('sbcisf', pvarout,'T',-1.) 813 814 CASE ( 'T' ) ! compute T in the top boundary layer at T- point 815 DO jj = 1,jpj 816 DO ji = 1,jpi 817 ikt = misfkt(ji,jj) 818 ikb = misfkb(ji,jj) 819 820 ! level fully include in the ice shelf boundary layer 821 DO jk = ikt, ikb - 1 822 ze3 = e3t_n(ji,jj,jk) 823 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 824 END DO 825 826 ! level partially include in ice shelf boundary layer 827 zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 828 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 829 END DO 830 END DO 831 END SELECT 832 ! 833 ! mask mean tbl value 834 pvarout(:,:) = pvarout(:,:) * ssmask(:,:) 835 ! 836 END SUBROUTINE sbc_isf_tbl 837 838 839 SUBROUTINE sbc_isf_div( phdivn ) 840 !!---------------------------------------------------------------------- 841 !! *** SUBROUTINE sbc_isf_div *** 842 !! 843 !! ** Purpose : update the horizontal divergence with the runoff inflow 844 !! 845 !! ** Method : 846 !! CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the 847 !! divergence and expressed in m/s 848 !! 849 !! ** Action : phdivn decreased by the runoff inflow 850 !!---------------------------------------------------------------------- 851 REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: phdivn ! horizontal divergence 852 ! 853 INTEGER :: ji, jj, jk ! dummy loop indices 854 INTEGER :: ikt, ikb 855 REAL(wp) :: zhk 856 REAL(wp) :: zfact ! local scalar 857 !!---------------------------------------------------------------------- 858 ! 859 zfact = 0.5_wp 860 ! 861 IF(.NOT.ln_linssh ) THEN ! need to re compute level distribution of isf fresh water 862 DO jj = 1,jpj 863 DO ji = 1,jpi 864 ikt = misfkt(ji,jj) 865 ikb = misfkt(ji,jj) 866 ! thickness of boundary layer at least the top level thickness 867 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt)) 868 869 ! determine the deepest level influenced by the boundary layer 870 DO jk = ikt, mbkt(ji,jj) 871 IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 872 END DO 873 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 874 misfkb(ji,jj) = ikb ! last wet level of the tbl 875 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 876 877 zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 878 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 879 END DO 880 END DO 881 END IF 882 ! 883 !== ice shelf melting distributed over several levels ==! 884 DO jj = 1,jpj 885 DO ji = 1,jpi 886 ikt = misfkt(ji,jj) 887 ikb = misfkb(ji,jj) 888 ! level fully include in the ice shelf boundary layer 889 DO jk = ikt, ikb - 1 890 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 891 & * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 892 END DO 893 ! level partially include in ice shelf boundary layer 894 phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 895 & + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 896 END DO 897 END DO 898 ! 899 END SUBROUTINE sbc_isf_div 900 179 ! 180 IF ( ln_isfcav_mlt ) THEN 181 ! 182 ! initialisation of cav variable 183 CALL isf_cav_init() 184 ! 185 ! read cav variable from restart 186 IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b) 187 ! 188 END IF 189 ! 190 IF ( ln_isfpar_mlt ) THEN 191 ! 192 ! initialisation of par variable 193 CALL isf_par_init() 194 ! 195 ! read par variable from restart 196 IF ( ln_rstart ) CALL isfrst_read('par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b) 197 ! 198 END IF 199 ! 200 END SUBROUTINE isf_mlt_init 901 201 !!====================================================================== 902 END MODULE sbcisf202 END MODULE isfmlt -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/LDF/ldfslp.F90
r10425 r11395 21 21 !!---------------------------------------------------------------------- 22 22 USE oce ! ocean dynamics and tracers 23 USE isf ! ice shelf 23 24 USE dom_oce ! ocean space and time domain 24 25 ! USE ldfdyn ! lateral diffusion: eddy viscosity coef. -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90
r10882 r11395 123 123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fmmflx !: freshwater budget: freezing/melting [Kg/m2/s] 124 124 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnf , rnf_b !: river runoff [Kg/m2/s] 125 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf , fwfisf_b !: ice shelf melting [Kg/m2/s]126 125 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwficb , fwficb_b !: iceberg melting [Kg/m2/s] 127 126 … … 174 173 & sfx (jpi,jpj) , sfx_b(jpi,jpj) , emp_tot(jpi,jpj), fmmflx(jpi,jpj), STAT=ierr(2) ) 175 174 ! 176 ALLOCATE( fwfisf (jpi,jpj),rnf (jpi,jpj) , sbc_tsc (jpi,jpj,jpts) , qsr_hc (jpi,jpj,jpk) , &177 & fwfisf_b(jpi,jpj),rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , &175 ALLOCATE( rnf (jpi,jpj) , sbc_tsc (jpi,jpj,jpts) , qsr_hc (jpi,jpj,jpk) , & 176 & rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , & 178 177 & fwficb (jpi,jpj), fwficb_b(jpi,jpj), STAT=ierr(3) ) 179 178 ! -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90
r10617 r11395 36 36 USE eosbn2 ! 37 37 USE sbcrnf , ONLY : l_rnfcpl 38 USE sbcisf , ONLY : l_isfcpl38 USE isf , ONLY : l_isfcpl, fwfisf_cav, fwfisf_par 39 39 #if defined key_cice 40 40 USE ice_domain_size, only: ncat … … 478 478 IF(lwp) WRITE(numout,*) 479 479 IF(lwp) WRITE(numout,*) ' iceshelf received from oasis ' 480 CALL ctl_stop('STOP','not coded') 480 481 ENDIF 481 482 ! … … 1405 1406 rnf(:,:) = rnf(:,:) + fwficb(:,:) ! iceberg added to runfofs 1406 1407 ENDIF 1407 IF( srcv(jpr_isf)%laction ) fwfisf (:,:) = - frcv(jpr_isf)%z3(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting)1408 IF( srcv(jpr_isf)%laction ) fwfisf_par(:,:) = - frcv(jpr_isf)%z3(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) 1408 1409 1409 1410 IF( ln_mixcpl ) THEN ; emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:) … … 1708 1709 ENDIF 1709 1710 IF( srcv(jpr_isf)%laction ) THEN ! iceshelf (fwfisf <0 mean melting) 1710 fwfisf (:,:) = - frcv(jpr_isf)%z3(:,:,1)1711 fwfisf_par(:,:) = - frcv(jpr_isf)%z3(:,:,1) 1711 1712 ENDIF 1712 1713 … … 1747 1748 ENDIF 1748 1749 IF( srcv(jpr_isf)%laction ) THEN ! iceshelf (fwfisf <0 mean melting) 1749 fwfisf (:,:) = - frcv(jpr_isf)%z3(:,:,1)1750 fwfisf_par(:,:) = - frcv(jpr_isf)%z3(:,:,1) 1750 1751 ENDIF 1751 1752 ! -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcfwb.F90
r10570 r11395 20 20 USE phycst ! physical constants 21 21 USE sbcrnf ! ocean runoffs 22 USE sbcisf! ice shelf melting contribution22 USE isf ! ice shelf melting contribution 23 23 USE sbcssr ! Sea-Surface damping terms 24 24 ! … … 104 104 ! 105 105 IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 106 y_fwfnow(1) = local_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf (:,:) - snwice_fmass(:,:) ) )106 y_fwfnow(1) = local_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) - snwice_fmass(:,:) ) ) 107 107 CALL mpp_delay_sum( 'sbcfwb', 'fwb', y_fwfnow(:), z_fwfprv(:), kt == nitend - nn_fsbc + 1 ) 108 108 z_fwfprv(1) = z_fwfprv(1) / area … … 159 159 ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:) 160 160 ! ! fwf global mean (excluding ocean to ice/snow exchanges) 161 z_fwf = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf (:,:) - snwice_fmass(:,:) ) ) / area161 z_fwf = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) - snwice_fmass(:,:) ) ) / area 162 162 ! 163 163 IF( z_fwf < 0._wp ) THEN ! spread out over >0 erp area to increase evaporation -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcmod.F90
r10499 r11395 37 37 #endif 38 38 USE sbcice_cice ! surface boundary condition: CICE sea-ice model 39 USE sbcisf ! surface boundary condition: ice-shelf40 39 USE sbccpl ! surface boundary condition: coupled formulation 41 40 USE cpl_oasis3 ! OASIS routines for coupling … … 43 42 USE sbcrnf ! surface boundary condition: runoffs 44 43 USE sbcapr ! surface boundary condition: atmo pressure 45 USE sbcisf! surface boundary condition: ice shelf44 USE isf ! surface boundary condition: ice shelf 46 45 USE sbcfwb ! surface boundary condition: freshwater budget 47 46 USE icbstp ! Icebergs … … 239 238 #endif 240 239 ! 241 IF( .NOT.ln_isf ) THEN !* No ice-shelf in the domain : allocate and set to zero242 IF( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )243 fwfisf (:,:) = 0._wp ; risf_tsc (:,:,:) = 0._wp244 fwfisf_b(:,:) = 0._wp ; risf_tsc_b(:,:,:) = 0._wp245 END IF246 240 IF( nn_ice == 0 ) THEN !* No sea-ice in the domain : ice fraction is always zero 247 241 IF( nn_components /= jp_iam_opa ) fr_i(:,:) = 0._wp ! except for OPA in SAS-OPA coupled case … … 329 323 IF( ln_ssr ) CALL sbc_ssr_init ! Sea-Surface Restoring initialization 330 324 ! 331 IF( ln_isf ) CALL sbc_isf_init ! Compute iceshelves332 !333 325 CALL sbc_rnf_init ! Runof initialization 334 326 ! … … 397 389 rnf_b (:,: ) = rnf (:,: ) 398 390 rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) 399 ENDIF400 IF( ln_isf ) THEN401 fwfisf_b (:,: ) = fwfisf (:,: )402 risf_tsc_b(:,:,:) = risf_tsc(:,:,:)403 391 ENDIF 404 392 ! … … 450 438 IF( .NOT. ln_passive_mode ) CALL lbc_lnk( 'sbcmod', emp, 'T', 1. ) ! ensure restartability with icebergs 451 439 ENDIF 452 453 IF( ln_isf ) CALL sbc_isf( kt ) ! compute iceshelves454 440 455 441 IF( ln_rnf ) CALL sbc_rnf( kt ) ! add runoffs to fresh water fluxes … … 554 540 ! 555 541 IF(ln_ctl) THEN ! print mean trends (used for debugging) 556 CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i- : ', mask1=tmask )557 CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf- : ', mask1=tmask )558 CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf- : ', mask1=tmask )542 CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask ) 543 CALL prt_ctl(tab2d_1=(emp-rnf) , clinfo1=' emp-rnf - : ', mask1=tmask ) 544 CALL prt_ctl(tab2d_1=(sfx-rnf) , clinfo1=' sfx-rnf - : ', mask1=tmask ) 559 545 CALL prt_ctl(tab2d_1=qns , clinfo1=' qns - : ', mask1=tmask ) 560 546 CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask ) -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcrnf.F90
r10523 r11395 19 19 USE phycst ! physical constants 20 20 USE sbc_oce ! surface boundary condition variables 21 USE sbcisf ! PM we could remove it I think21 USE isf ! ice shelf 22 22 USE eosbn2 ! Equation Of State 23 23 USE closea ! closed seas -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90
r10425 r11395 28 28 USE sbc_oce ! surface boundary condition: ocean 29 29 USE sbcrnf ! river runoffs 30 USE sbcisf! ice shelf melting30 USE isf ! ice shelf melting 31 31 USE zdf_oce ! ocean vertical mixing 32 32 USE domvvl ! variable volume … … 312 312 ztc_f = ztc_n + atfp * ztc_d 313 313 ! 314 IF( jk == mikt(ji,jj) ) THEN ! first level 315 ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj) - emp(ji,jj) ) & 316 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) 314 IF( jk == 1 ) THEN ! first level 315 ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) ) 317 316 ztc_f = ztc_f - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 318 317 ENDIF 319 318 IF( ln_rnf_depth ) THEN 320 319 ! Rivers are not just at the surface must go down to nk_rnf(ji,jj) 321 IF( mikt(ji,jj) <=jk .and.jk <= nk_rnf(ji,jj) ) THEN320 IF( jk <= nk_rnf(ji,jj) ) THEN 322 321 ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj) - rnf(ji,jj) ) ) & 323 322 & * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) 324 323 ENDIF 325 324 ELSE 326 IF( jk == mikt(ji,jj)) THEN ! first level325 IF( jk == 1 ) THEN ! first level 327 326 ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj) - rnf(ji,jj) ) ) 328 327 ENDIF … … 341 340 ! ice shelf 342 341 IF( ll_isf ) THEN 343 ! level fully include in the Losch_2008 ice shelf boundary layer 344 IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) ) & 345 ztc_f = ztc_f - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) ) & 346 & * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) 347 ! level partially include in Losch_2008 ice shelf boundary layer 348 IF ( jk == misfkb(ji,jj) ) & 349 ztc_f = ztc_f - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) ) & 350 & * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj) 342 IF ( ln_isfcav_mlt ) THEN 343 ! level fully include in the Losch_2008 ice shelf boundary layer 344 IF ( jk >= misfkt_cav(ji,jj) .AND. jk < misfkb_cav(ji,jj) ) THEN 345 ztc_f = ztc_f - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) ) & 346 & * e3t_n(ji,jj,jk) / rhisf_tbl_cav(ji,jj) 347 ze3t_f = ze3t_f - zfact2 * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) ) & 348 & * e3t_n(ji,jj,jk) / rhisf_tbl_cav(ji,jj) 349 END IF 350 ! level partially include in Losch_2008 ice shelf boundary layer 351 IF ( jk == misfkb_cav(ji,jj) ) THEN 352 ztc_f = ztc_f - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) ) & 353 & * e3t_n(ji,jj,jk) / rhisf_tbl_cav(ji,jj) * rfrac_tbl_cav(ji,jj) 354 ze3t_f = ze3t_f - zfact2 * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) ) & 355 & * e3t_n(ji,jj,jk) / rhisf_tbl_cav(ji,jj) * rfrac_tbl_cav(ji,jj) 356 END IF 357 END IF 358 IF ( ln_isfpar_mlt ) THEN 359 ! level fully include in the Losch_2008 ice shelf boundary layer 360 IF ( jk >= misfkt_par(ji,jj) .AND. jk < misfkb_par(ji,jj) ) THEN 361 ztc_f = ztc_f - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) ) & 362 & * e3t_n(ji,jj,jk) / rhisf_tbl_par(ji,jj) 363 ze3t_f = ze3t_f - zfact2 * ( fwfisf_par_b(ji,jj) - fwfisf_par(ji,jj) ) & 364 & * e3t_n(ji,jj,jk) / rhisf_tbl_par(ji,jj) 365 END IF 366 ! level partially include in Losch_2008 ice shelf boundary layer 367 IF ( jk == misfkb_par(ji,jj) ) THEN 368 ztc_f = ztc_f - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) ) & 369 & * e3t_n(ji,jj,jk) / rhisf_tbl_par(ji,jj) * rfrac_tbl_par(ji,jj) 370 ze3t_f = ze3t_f - zfact2 * ( fwfisf_par_b(ji,jj) - fwfisf_par(ji,jj) ) & 371 & * e3t_n(ji,jj,jk) / rhisf_tbl_par(ji,jj) * rfrac_tbl_par(ji,jj) 372 END IF 373 END IF 351 374 END IF 352 375 ! -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/trasbc.F90
r10499 r11395 22 22 USE sbcmod ! ln_rnf 23 23 USE sbcrnf ! River runoff 24 USE sbcisf ! Ice shelf25 USE iscplini ! Ice sheet coupling26 24 USE traqsr ! solar radiation penetration 27 25 USE trd_oce ! trends: ocean variables … … 155 153 ! 156 154 !---------------------------------------- 157 ! Ice Shelf effects (ISF)158 ! tbl treated as in Losh (2008) JGR159 !----------------------------------------160 !161 !!gm BUG ? Why no differences between non-linear and linear free surface ?162 !!gm probably taken into account in r1_hisf_tbl : to be verified163 IF( ln_isf ) THEN164 zfact = 0.5_wp165 DO jj = 2, jpj166 DO ji = fs_2, fs_jpim1167 !168 ikt = misfkt(ji,jj)169 ikb = misfkb(ji,jj)170 !171 ! level fully include in the ice shelf boundary layer172 ! sign - because fwf sign of evapo (rnf sign of precip)173 DO jk = ikt, ikb - 1174 ! compute trend175 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) &176 & + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) ) &177 & * r1_hisf_tbl(ji,jj)178 END DO179 180 ! level partially include in ice shelf boundary layer181 ! compute trend182 tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem) &183 & + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) ) &184 & * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)185 186 END DO187 END DO188 END IF189 !190 !----------------------------------------191 155 ! River Runoff effects 192 156 !---------------------------------------- … … 242 206 #endif 243 207 ! 244 !----------------------------------------245 ! Ice Sheet coupling imbalance correction to have conservation246 !----------------------------------------247 !248 IF( ln_iscpl .AND. ln_hsb) THEN ! input of heat and salt due to river runoff249 DO jk = 1,jpk250 DO jj = 2, jpj251 DO ji = fs_2, fs_jpim1252 zdep = 1._wp / e3t_n(ji,jj,jk)253 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep254 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep255 END DO256 END DO257 END DO258 ENDIF259 260 208 IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics 261 209 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ZDF/zdfmxl.F90
r10425 r11395 12 12 !!---------------------------------------------------------------------- 13 13 USE oce ! ocean dynamics and tracers variables 14 USE isf ! ice shelf 14 15 USE dom_oce ! ocean space and time domain variables 15 16 USE trc_oce , ONLY: l_offline ! ocean space and time domain variables -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90
r10588 r11395 60 60 USE diacfl ! CFL diagnostics (dia_cfl_init routine) 61 61 USE step ! NEMO time-stepping (stp routine) 62 USE isfmlt ! ice shelf (isf_mlt_init routine) 62 63 USE icbini ! handle bergs, initialisation 63 64 USE icbstp ! handle bergs, calving, themodynamics and transport … … 452 453 ! ! Active tracers 453 454 IF( ln_traqsr ) CALL tra_qsr_init ! penetrative solar radiation qsr 455 IF (ln_isf ) CALL isf_mlt_init ! ice shelf 454 456 CALL tra_bbc_init ! bottom heat flux 455 457 CALL tra_bbl_init ! advective (and/or diffusive) bottom boundary layer scheme -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90
r11287 r11395 113 113 IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 114 114 IF( ln_bdy ) CALL bdy_dta ( kstp, time_offset=+1 ) ! update dynamic & tracer data at open boundaries 115 IF( ln_isf ) CALL isf_mlt ( kstp ) 115 116 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 116 117 … … 240 241 CALL tra_sbc ( kstp ) ! surface boundary condition 241 242 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 243 IF( ln_isf ) CALL tra_isf ( kstp ) ! ice shelf heat flux 242 244 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 243 245 IF( ln_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme -
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step_oce.F90
r10068 r11395 22 22 USE sbcwave ! Wave intialisation 23 23 24 USE isfmlt ! ice shelf boundary condition (isf_mlt routine) 25 24 26 USE traqsr ! solar radiation penetration (tra_qsr routine) 27 USE traisf ! ice shelf (tra_isf routine) 25 28 USE trasbc ! surface boundary condition (tra_sbc routine) 26 29 USE trabbc ! bottom boundary condition (tra_bbc routine)
Note: See TracChangeset
for help on using the changeset viewer.