Changeset 8498 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM
- Timestamp:
- 2017-09-05T19:53:41+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r8486 r8498 332 332 333 333 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: afx_tot !: ice concentration tendency (total) [s-1] 334 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: afx_thd !: ice concentration tendency (thermodynamics) [s-1]335 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: afx_dyn !: ice concentration tendency (dynamics) [s-1]336 334 337 335 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_bog !: salt flux due to ice growth/melt [PSU/m2/s] … … 479 477 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_dmi_dyn !: Change in ice mass due to ice dynamics (kg/m2/s) 480 478 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_dms_dyn !: Change in snow mass due to ice dynamics (kg/m2/s) 481 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_xmtrp_ice !: X-component of ice mass transport (kg/s)482 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_ymtrp_ice !: Y-component of ice mass transport (kg/s)483 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_xmtrp_snw !: X-component of snow mass transport (kg/s)484 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_ymtrp_snw !: Y-component of snow mass transport (kg/s)485 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_xatrp !: X-component of area transport (m2/s)486 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_yatrp !: Y-component of area transport (m2/s)487 479 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_fc_bo !: Bottom conduction flux (W/m2) 488 480 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_fc_su !: Surface conduction flux (W/m2) 489 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_utau_oi !: X-direction ocean-ice stress490 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vtau_oi !: Y-direction ocean-ice stress491 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_dssh_dx !: X-direction sea-surface tilt term (N/m2)492 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_dssh_dy !: X-direction sea-surface tilt term (N/m2)493 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_corstrx !: X-direction coriolis stress (N/m2)494 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_corstry !: Y-direction coriolis stress (N/m2)495 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_intstrx !: X-direction internal stress (N/m2)496 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_intstry !: Y-direction internal stress (N/m2)497 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_sig1 !: Average normal stress in sea ice498 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_sig2 !: Maximum shear stress in sea ice499 481 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_shear !: Maximum shear of sea-ice velocity field 500 482 … … 534 516 & wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) , & 535 517 & wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) , & 536 & afx_tot(jpi,jpj) , afx_thd(jpi,jpj), afx_dyn(jpi,jpj) , rn_amax_2d(jpi,jpj),&518 & afx_tot(jpi,jpj) , rn_amax_2d(jpi,jpj), & 537 519 & fhtur (jpi,jpj) , qlead (jpi,jpj) , & 538 520 & sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) , & … … 616 598 ALLOCATE( t_si (jpi,jpj,jpl) , tm_si(jpi,jpj) , & 617 599 diag_dmi_dyn(jpi,jpj) , diag_dms_dyn(jpi,jpj) , & 618 diag_xmtrp_ice(jpi,jpj), diag_ymtrp_ice(jpi,jpj), &619 diag_xmtrp_snw(jpi,jpj), diag_ymtrp_snw(jpi,jpj), &620 diag_xatrp(jpi,jpj) , diag_yatrp(jpi,jpj) , &621 600 diag_fc_bo(jpi,jpj) , diag_fc_su(jpi,jpj) , & 622 diag_utau_oi(jpi,jpj) , diag_vtau_oi(jpi,jpj) , &623 diag_dssh_dx(jpi,jpj) , diag_dssh_dy(jpi,jpj) , &624 diag_corstrx(jpi,jpj) , diag_corstry(jpi,jpj) , &625 diag_intstrx(jpi,jpj) , diag_intstry(jpi,jpj) , &626 diag_sig1(jpi,jpj) , diag_sig2(jpi,jpj) , &627 601 diag_shear(jpi,jpj) , STAT = ierr(ii) ) 628 602 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90
r8486 r8498 27 27 REAL(wp), PUBLIC, PARAMETER :: rn_alb_oce = 0.066 !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 28 28 29 !!gm real variable name starting with a "c" NOT DOCTOR !!!!30 29 INTEGER :: albd_init = 0 ! control flag for initialization 31 REAL(wp) , PARAMETER :: c1= 0.05 ! snow thickness (only for nn_ice_alb=0)32 REAL(wp) , PARAMETER :: c2= 0.10 ! " "33 REAL(wp) , PARAMETER :: rcloud 34 REAL(wp) , PARAMETER :: r1_c1 = 1. / c135 REAL(wp) , PARAMETER :: r1_c2 = 1. / c230 REAL(wp) , PARAMETER :: rc1 = 0.05 ! snow thickness (only for nn_ice_alb=0) 31 REAL(wp) , PARAMETER :: rc2 = 0.10 ! " " 32 REAL(wp) , PARAMETER :: rcloud = 0.06 ! cloud effect on albedo (only-for nn_ice_alb=0) 33 REAL(wp) , PARAMETER :: r1_c1 = 1. / rc1 34 REAL(wp) , PARAMETER :: r1_c2 = 1. / rc2 36 35 37 36 ! !!* namelist namsbc_alb * … … 114 113 IF ( .NOT. ld_pnd ) THEN !--- No melt ponds OR radiatively inactive melt ponds 115 114 ! Bare ice albedo is prescribed, with implicit assumption on pond fraction and depth 116 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice ) ; zalb(:,:,:) = rn_alb_imlt 117 ! !!! MV I think we could replace rt0_ice by rt0 and get rid of rt0 118 ELSE WHERE ; zalb(:,:,:) = rn_alb_idry 119 END WHERE 115 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 ) ; zalb(:,:,:) = rn_alb_imlt 116 ELSEWHERE ; zalb(:,:,:) = rn_alb_idry 117 END WHERE 120 118 ENDIF 121 119 … … 127 125 ! 128 126 ! ! Thickness-dependent bare ice albedo 129 WHERE ( 1.5 < ph_ice ) ;zalb_it = zalb130 ELSE WHERE( 1.0 < ph_ice .AND. ph_ice <= 1.5 ) ;zalb_it = 0.472 + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )131 ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 ) ;zalb_it = 0.2467 + 0.7049 * ph_ice &132 & - 0.8608 * ph_ice * ph_ice &133 & + 0.3812 * ph_ice * ph_ice * ph_ice134 ELSE WHERE ;zalb_it = 0.1 + 3.6 * ph_ice127 WHERE ( 1.5 < ph_ice ) ; zalb_it = zalb 128 ELSEWHERE( 1.0 < ph_ice .AND. ph_ice <= 1.5 ) ; zalb_it = 0.472 + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) 129 ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 ) ; zalb_it = 0.2467 + 0.7049 * ph_ice & 130 & - 0.8608 * ph_ice * ph_ice & 131 & + 0.3812 * ph_ice * ph_ice * ph_ice 132 ELSEWHERE ; zalb_it = 0.1 + 3.6 * ph_ice 135 133 END WHERE 136 134 ! … … 139 137 zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd ) 140 138 ! ! Snow-free ice albedo is a function of pond fraction 141 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 _ice)139 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 ) 142 140 zalb_it = zalb_it * ( 1. - pafrac_pnd ) + zalb_pnd * pafrac_pnd 143 141 END WHERE … … 150 148 DO ji = 1, jpi 151 149 ! ! Freezing snow 152 ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2153 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )150 ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > rc2 151 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - rc1 ) ) ) 154 152 zalb_sf = ( 1._wp - zswitch ) * ( zalb_it(ji,jj,jl) & 155 153 & + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * r1_c1 ) & … … 157 155 ! 158 156 ! ! Melting snow 159 ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2160 zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )157 ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > rc2 158 zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - rc2 ) ) 161 159 zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * r1_c2 ) & 162 160 & + zswitch * rn_alb_smlt … … 211 209 ! 212 210 ! ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions 213 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 _ice)211 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 ) 214 212 zalb_it = zalb_it * ( 1. - pafrac_pnd ) + zalb_pnd * pafrac_pnd 215 213 END WHERE … … 256 254 ! !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time) 257 255 zafrac_pnd = 0._wp 258 IF ( ld_pnd ) THEN 259 WHERE( ph_snw == 0._wp ) ; zafrac_pnd = pafrac_pnd ; END WHERE! Snow fully "shades" melt ponds256 IF ( ld_pnd ) THEN 257 WHERE( ph_snw == 0._wp ) zafrac_pnd = pafrac_pnd ! Snow fully "shades" melt ponds 260 258 ENDIF 261 259 ! 262 260 ! !--- Specific snow fraction (for now, prescribed) 263 WHERE ( ph_snw > 0._wp ) ;zafrac_snw = 1.264 ELSE WHERE ;zafrac_snw = 0.261 WHERE( ph_snw > 0._wp ) ; zafrac_snw = 1. 262 ELSEWHERE ; zafrac_snw = 0. 265 263 END WHERE 266 264 ! … … 273 271 z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 274 272 z1_c2 = 1. / 0.05 275 WHERE ( 1.5 < ph_ice ) ;zalb_ice = zalb276 ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 ) ; zalb_ice = zalb + ( 0.18 - zalb) * z1_c1 * &273 WHERE ( 1.5 < ph_ice ) ; zalb_ice = zalb 274 ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 ) ; zalb_ice = zalb + ( 0.18 - zalb ) * z1_c1 * & 277 275 & ( LOG(1.5) - LOG(ph_ice) ) 278 ELSE WHERE ;zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice276 ELSEWHERE ; zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice 279 277 END WHERE 280 278 ! … … 282 280 z1_c2 = 1. / 0.03 283 281 ! 284 WHERE( pt_ice < rt0_snow ) ; zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 );285 ELSE WHERE ; zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 );282 WHERE( pt_ice < rt0_snow ) ; zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 ) 283 ELSEWHERE ; zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 ) 286 284 END WHERE 287 285 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90
r8491 r8498 26 26 USE lib_mpp ! MPP library 27 27 USE timing ! Timing 28 USE iom ! 28 29 29 30 IMPLICIT NONE … … 53 54 INTEGER :: ji, jj, jk, jl ! dummy loop indices 54 55 REAL(wp) :: zsal, zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b, zzc 56 REAL(wp), DIMENSION(jpi,jpj) :: zafx ! concentration trends diag 55 57 !!---------------------------------------------------------------------- 56 58 IF( nn_timing == 1 ) CALL timing_start('icecor') … … 65 67 ! 66 68 ! 67 ! !----------------------------------------------------- 68 IF( kn == 2 ) THEN ! thickness of the smallest category above himin ! 69 IF( kn == 2 ) THEN 69 70 ! !----------------------------------------------------- 70 ! 71 DO jj = 1, jpj 72 DO ji = 1, jpi 73 !!gm replace this 74 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) ) !0 if no ice and 1 if yes 75 ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch 76 !!gm by more readable coding (not slower coding since already a IF in the loop): 77 ! IF( a_i(ji,jj,1) >= epsi20 ) ht_i(ji,jj,1) = v_i (ji,jj,1) / a_i(ji,jj,1) 78 !!gm 79 IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN 80 a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin 81 ENDIF 82 END DO 83 END DO 71 ! ! thickness of the smallest category above himin ! 72 ! !----------------------------------------------------- 73 WHERE( a_i(:,:,1) >= epsi20 ) ; ht_i(:,:,1) = v_i (:,:,1) / a_i(:,:,1) 74 ELSEWHERE ; ht_i(:,:,1) = 0._wp 75 END WHERE 76 WHERE( ht_i(:,:,1) < rn_himin ) a_i (:,:,1) = a_i (:,:,1) * ht_i(:,:,1) / rn_himin 84 77 ! 85 78 ENDIF 86 87 79 ! !----------------------------------------------------- 88 at_i(:,:) = a_i(:,:,1) ! ice concentration should not exceed amax ! 89 DO jl = 2, jpl !----------------------------------------------------- 90 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 80 ! ! ice concentration should not exceed amax ! 81 ! !----------------------------------------------------- 82 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 83 DO jl = 1, jpl 84 WHERE( at_i(:,:) > rn_amax_2d(:,:) ) a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:) 91 85 END DO 92 !93 !!gm Question it seams to me that we have the following equality (dropping the "(ji,jj)":94 ! ( 1. - ( 1. - rn_amax_2d / at_i ) ) = ( 1. - ( at_i - rn_amax_2d ) / at_i )95 ! = ( at_i - ( at_i - rn_amax_2d ) ) / at_i96 ! = ( + rn_amax_2d ) / at_i97 ! = rn_amax_2d / at_i98 ! No ? if yes see "!!gm better" juste below99 !gm100 DO jl = 1, jpl101 DO jj = 1, jpj102 DO ji = 1, jpi103 IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN104 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )105 !!gm better: a_i(ji,jj,jl) = a_i(ji,jj,jl) * rn_amax_2d(ji,jj) / at_i(ji,jj)106 ENDIF107 END DO108 END DO109 END DO110 !!gm Other question: why testing a_i(ji,jj,jl) > 0._wp ? a_i is >=0, a multiplication by 0 does not change the results....111 !!gm so at the end, the loop can be recoded without IF as:112 ! WHERE( at_i(:,:) > rn_amax_2d(:,:) )113 ! DO jl = 1, jpl114 ! a_i(:,:,jl) = a_i(:,:,jl) * MAX( rn_amax_2d(:,:), at_i(:,:) ) / at_i(:,:)115 ! END DO116 ! END WHERE117 !!gm No?118 86 119 87 ! !----------------------------------------------------- 120 IF ( nn_icesal == 2 ) THEN ! Ice salinity bounds!88 IF ( nn_icesal == 2 ) THEN ! salinity stays in bounds [Simin,Simax] ! 121 89 ! !----------------------------------------------------- 122 90 zzc = rhoic * r1_rdtice 123 DO jl = 1, jpl ! salinity stays in bounds [Simin,Simax]91 DO jl = 1, jpl 124 92 DO jj = 1, jpj 125 93 DO ji = 1, jpi 126 IF( v_i(ji,jj,jl) > 0._wp ) THEN ! clem: useless IF ??? 127 zsal = smv_i(ji,jj,jl) 128 smv_i(ji,jj,jl) = MIN( MAX( rn_simin*v_i(ji,jj,jl) , smv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl) ) 129 ! associated salt flux 130 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * zzc 131 ENDIF 94 zsal = smv_i(ji,jj,jl) 95 smv_i(ji,jj,jl) = MIN( MAX( rn_simin*v_i(ji,jj,jl) , smv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl) ) 96 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * zzc ! associated salt flux 132 97 END DO 133 98 END DO … … 157 122 END DO 158 123 CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. ) ! lateral boundary conditions 159 !160 !!gm I think masking here is unnecessary, u_ice already masked and we only introduce zeros in the field161 u_ice(:,:) = u_ice(:,:) * umask(:,:,1) ! mask velocities162 v_ice(:,:) = v_ice(:,:) * vmask(:,:,1)163 124 ENDIF 164 125 … … 172 133 ! !----------------------------------------------------- 173 134 CASE( 1 ) !--- dyn trend diagnostics 174 DO jl = 1, jpl175 afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice176 END DO177 135 ! 178 136 !!gm here I think the number of ice cat is too small to use a SUM instruction... … … 188 146 END DO 189 147 END DO 148 ! ! concentration tendency (dynamics) 149 zafx (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 150 afx_tot(:,:) = zafx(:,:) 151 IF( iom_use('afxdyn') ) CALL iom_put( 'afxdyn' , zafx(:,:) ) 190 152 ! 191 153 CASE( 2 ) !--- thermo trend diagnostics & ice aging 192 154 ! 193 DO jl = 1, jpl 194 oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice ! ice natural aging incrementation 195 afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice ! thermo tendancy 196 END DO 197 afx_tot(:,:) = afx_thd(:,:) + afx_dyn(:,:) 155 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice ! ice natural aging incrementation 198 156 ! 199 157 !!gm here I think the number of ice cat is too small to use a SUM instruction... … … 209 167 END DO 210 168 END DO 169 ! ! concentration tendency (total + thermo) 170 zafx (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 171 afx_tot(:,:) = afx_tot(:,:) + zafx(:,:) 172 IF( iom_use('afxthd') ) CALL iom_put( 'afxthd' , zafx(:,:) ) 173 IF( iom_use('afxtot') ) CALL iom_put( 'afxtot' , afx_tot(:,:) ) 211 174 ! 212 175 END SELECT -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
r8486 r8498 71 71 REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt 72 72 !!---------------------------------------------------------------------- 73 74 !!gm Note that glo_sum for a 2D or 3D array use a multiplication by tmask_i(ji,jj) 75 !! so below the * tmask(:,:,1) is useless ===>> I have removed them 76 !! I also move the conversion factor from then glo_sum argument (become a single multiplication 77 73 ! 78 74 IF( icount == 0 ) THEN 79 75 ! ! salt flux … … 271 267 !WRITE(numout,*) ' sea-ice stress utau_ice : ', utau_ice(ji,jj) 272 268 !WRITE(numout,*) ' sea-ice stress vtau_ice : ', vtau_ice(ji,jj) 273 !WRITE(numout,*) ' oceanic speed u : ', u_oce(ji,jj)274 !WRITE(numout,*) ' oceanic speed v : ', v_oce(ji,jj)275 269 !WRITE(numout,*) ' sst : ', sst_m(ji,jj) 276 270 !WRITE(numout,*) ' sss : ', sss_m(ji,jj) … … 556 550 WRITE(numout,*) ' utau : ', utau (ji,jj) 557 551 WRITE(numout,*) ' vtau : ', vtau (ji,jj) 558 WRITE(numout,*) ' oc. vel. u : ', u_oce (ji,jj)559 WRITE(numout,*) ' oc. vel. v : ', v_oce (ji,jj)560 552 ENDIF 561 553 … … 669 661 CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', tab2d_2=vtau , clinfo2= ' vtau : ') 670 662 CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ') 671 CALL prt_ctl(tab2d_1=u_oce , clinfo1= ' u_oce : ', tab2d_2=v_oce , clinfo2= ' v_oce : ')672 663 673 664 END SUBROUTINE ice_prt3D -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
r8486 r8498 65 65 IF( kt == nit000 .AND. lwp ) THEN 66 66 WRITE(numout,*) 67 WRITE(numout,*)'ice_forcing_tau '67 WRITE(numout,*)'ice_forcing_tau : Surface boundary condition for sea ice (momentum)' 68 68 WRITE(numout,*)'~~~~~~~~~~~~~~~' 69 69 ENDIF … … 124 124 IF( kt == nit000 .AND. lwp ) THEN 125 125 WRITE(numout,*) 126 WRITE(numout,*)'ice_forcing_flx '126 WRITE(numout,*)'ice_forcing_flx : Surface boundary condition for sea ice (flux)' 127 127 WRITE(numout,*)'~~~~~~~~~~~~~~~' 128 128 ENDIF … … 132 132 133 133 ! albedo depends on cloud fraction because of non-linear spectral effects 134 DO jl = 1, jpl135 DO jj = 2, jpjm1136 DO ji = 2, jpim1137 134 !!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument ! 138 alb_ice(ji,jj,jl) = ( 1. - cldf_ice ) * zalb_cs(ji,jj,jl) + cldf_ice * zalb_os(ji,jj,jl) 139 END DO 140 END DO 141 END DO 142 CALL lbc_lnk( alb_ice, 'T', 1. ) 135 alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 143 136 144 137 SELECT CASE( ksbc ) !== fluxes over sea ice ==! … … 198 191 INTEGER :: jl ! dummy loop index 199 192 ! 200 REAL(wp), DIMENSION(jpi,jpj) :: zalb_m ! Mean albedo over all categories 201 REAL(wp), DIMENSION(jpi,jpj) :: ztem_m ! Mean temperature over all categories 202 ! 203 REAL(wp), DIMENSION(jpi,jpj) :: z_qsr_m ! Mean solar heat flux over all categories 204 REAL(wp), DIMENSION(jpi,jpj) :: z_qns_m ! Mean non solar heat flux over all categories 205 REAL(wp), DIMENSION(jpi,jpj) :: z_evap_m ! Mean sublimation over all categories 206 REAL(wp), DIMENSION(jpi,jpj) :: z_dqn_m ! Mean d(qns)/dT over all categories 207 REAL(wp), DIMENSION(jpi,jpj) :: z_devap_m ! Mean d(evap)/dT over all categories 193 REAL(wp), DIMENSION(jpi,jpj) :: z1_at_i ! inverse of concentration 194 ! 195 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_qsr_m ! Mean solar heat flux over all categories 196 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_qns_m ! Mean non solar heat flux over all categories 197 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_evap_m ! Mean sublimation over all categories 198 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_dqn_m ! Mean d(qns)/dT over all categories 199 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories 200 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zalb_m ! Mean albedo over all categories 201 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztem_m ! Mean temperature over all categories 208 202 !!---------------------------------------------------------------------- 209 203 ! 210 204 IF( nn_timing == 1 ) CALL timing_start('ice_flx_dist') 211 205 ! 206 WHERE ( at_i (:,:) > 0._wp ) ; z1_at_i(:,:) = 1._wp / at_i (:,:) 207 ELSEWHERE ; z1_at_i(:,:) = 0._wp 208 END WHERE 209 212 210 SELECT CASE( k_limflx ) !== averaged on all ice categories ==! 213 211 ! 214 212 CASE( 0 , 1 ) 215 z_qns_m (:,:) = fice_ice_ave( pqns_ice (:,:,:) ) 216 z_qsr_m (:,:) = fice_ice_ave( pqsr_ice (:,:,:) ) 217 z_dqn_m (:,:) = fice_ice_ave( pdqn_ice (:,:,:) ) 218 z_evap_m (:,:) = fice_ice_ave( pevap_ice (:,:,:) ) 219 z_devap_m(:,:) = fice_ice_ave( pdevap_ice(:,:,:) ) 220 !!gm faster coding 221 ! REAL(wp), DIMENSION(jpi,jpj) :: z1_at_i ! 222 ! ... 223 ! WHERE ( at_i (:,:) > 0._wp ) ; z1_at_i(:,:) = 1._wp / at_i (:,:) 224 ! ELSEWHERE ; z1_at_i(:,:) = 0._wp 225 ! END WHERE 226 ! z_qns_m (:,:) = SUM( a_i(:,:,:) * pqns_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 227 ! z_qsr_m (:,:) = SUM( a_i(:,:,:) * pqsr_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 228 ! z_dqn_m (:,:) = SUM( a_i(:,:,:) * pdqn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 229 ! z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 230 ! z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:) 231 !! and remove the 2 functions: fice_ice_ave and fice_cell_ave 232 !!gm 213 ! 214 ALLOCATE( z_qns_m(jpi,jpj), z_qsr_m(jpi,jpj), z_dqn_m(jpi,jpj), z_evap_m(jpi,jpj), z_devap_m(jpi,jpj) ) 215 ! 216 z_qns_m (:,:) = SUM( a_i(:,:,:) * pqns_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 217 z_qsr_m (:,:) = SUM( a_i(:,:,:) * pqsr_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 218 z_dqn_m (:,:) = SUM( a_i(:,:,:) * pdqn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 219 z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 220 z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:) 233 221 DO jl = 1, jpl 234 pdqn_ice (:,:,jl) = z_dqn_m (:,:)235 pdevap_ice(:,:,jl) = z_devap_m(:,:)236 222 pqns_ice (:,:,jl) = z_qns_m (:,:) 237 223 pqsr_ice (:,:,jl) = z_qsr_m (:,:) 224 pdqn_ice (:,:,jl) = z_dqn_m (:,:) 238 225 pevap_ice (:,:,jl) = z_evap_m(:,:) 226 pdevap_ice(:,:,jl) = z_devap_m(:,:) 239 227 END DO 240 228 ! 229 DEALLOCATE( z_qns_m, z_qsr_m, z_dqn_m, z_evap_m, z_devap_m ) 230 ! 241 231 END SELECT 242 232 ! 243 233 SELECT CASE( k_limflx ) !== redistribution on all ice categories ==! 234 ! 244 235 CASE( 1 , 2 ) 245 236 ! 246 zalb_m(:,:) = fice_ice_ave( palb_ice(:,:,:) ) 247 ztem_m(:,:) = fice_ice_ave( ptn_ice (:,:,:) ) 237 ALLOCATE( zalb_m(jpi,jpj), ztem_m(jpi,jpj) ) 238 ! 239 zalb_m(:,:) = SUM( a_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:) 240 ztem_m(:,:) = SUM( a_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 248 241 DO jl = 1, jpl 249 242 pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) ) … … 252 245 END DO 253 246 ! 247 DEALLOCATE( zalb_m, ztem_m ) 248 ! 254 249 END SELECT 255 250 ! … … 257 252 ! 258 253 END SUBROUTINE ice_flx_dist 259 260 !!gm TO BE REMOVED ====>>>>>261 FUNCTION fice_cell_ave ( ptab )262 !!--------------------------------------------------------------------------263 !! * Compute average over categories, for grid cell (ice covered and free ocean)264 !!--------------------------------------------------------------------------265 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) :: ptab266 REAL(wp), DIMENSION(jpi,jpj) :: fice_cell_ave267 INTEGER :: jl268 !!--------------------------------------------------------------------------269 fice_cell_ave(:,:) = a_i(:,:,1) * ptab (:,:,1)270 DO jl = 2, jpl271 fice_cell_ave(:,:) = fice_cell_ave(:,:) + a_i(:,:,jl) * ptab (:,:,jl)272 END DO273 END FUNCTION fice_cell_ave274 275 276 FUNCTION fice_ice_ave ( ptab )277 !!--------------------------------------------------------------------------278 !! * Compute average over categories, for ice covered part of grid cell279 !!--------------------------------------------------------------------------280 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) :: ptab !281 REAL(wp), DIMENSION(jpi,jpj) :: fice_ice_ave282 !!--------------------------------------------------------------------------283 WHERE ( at_i (:,:) > 0.0_wp ) ; fice_ice_ave (:,:) = fice_cell_ave( ptab (:,:,:) ) / at_i (:,:)284 ELSEWHERE ; fice_ice_ave (:,:) = 0.0_wp285 END WHERE286 END FUNCTION fice_ice_ave287 288 !!gm <<<<<<==== end of TO BE REMOVED289 254 290 255 #else -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90
r8486 r8498 546 546 WRITE(numout,*) 547 547 WRITE(numout,*) 'ice_ist_init : ice parameters inititialisation ' 548 WRITE(numout,*) '~~~~~~~~~~~~~~~' 549 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_limini = ', ln_limini 550 WRITE(numout,*) ' ice initialization from a netcdf file ln_limini_file = ', ln_limini_file 551 WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst 552 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n 553 WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s 554 WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n 555 WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s 556 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n 557 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s 558 WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n 559 WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s 560 WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n 561 WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s 548 WRITE(numout,*) '~~~~~~~~~~~~' 549 WRITE(numout,*) ' Namelist namiceini' 550 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_limini = ', ln_limini 551 WRITE(numout,*) ' ice initialization from a netcdf file ln_limini_file = ', ln_limini_file 552 WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst 553 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n 554 WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s 555 WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n 556 WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s 557 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n 558 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s 559 WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n 560 WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s 561 WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n 562 WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s 562 563 ENDIF 563 564 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
r8486 r8498 426 426 ENDIF 427 427 ! 428 !!gm use of rswitch not faster as there is already IF in the DO-loop ==>>> use IF which is more readable429 ! rswitch = MAX( 0._wp , SIGN( 1._wp , v_i_2d(ji,jl1) - epsi10 ) )430 ! zworkv(ji) = pdvice(ji,jl) / MAX( v_i_2d(ji,jl1), epsi10 ) * rswitch431 ! !432 ! rswitch = MAX( 0._wp , SIGN( 1._wp , a_i_2d(ji,jl1) - epsi10 ) )433 ! zworka(ji) = pdaice(ji,jl) / MAX( a_i_2d(ji,jl1), epsi10 ) * rswitch434 435 428 IF( v_i_2d(ji,jl1) >= epsi10 ) THEN ; zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1) 436 429 ELSE ; zworkv(ji) = 0._wp … … 439 432 ELSE ; zworka(ji) = 0._wp 440 433 ENDIF 441 !!gm end442 434 ! 443 435 a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) ! Ice areas -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
r8486 r8498 48 48 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging 49 49 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting 50 50 ! 51 51 REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier 52 52 REAL(wp), PARAMETER :: kraft = 0.5_wp ! rafting multipliyer 53 54 !!gm Cp is 1) not DOCTOR, 55 !! 2) misleading name: heat capacity instead of a constant, 56 !! 3) recomputed at each time-step, whereas it is stored in the module memory... 57 !! ===>>> compute it one for all inside the IF( kt == nit000 ) (i.e. without the ".AND. lwp") 58 REAL(wp) :: Cp ! ??? !!gm Not doctor ! 59 53 ! 54 REAL(wp) :: zdrho ! 60 55 ! 61 56 ! … … 111 106 INTEGER :: niter ! local integer 112 107 INTEGER :: iterate_ridging ! if =1, repeat the ridging 113 REAL(wp) :: za, zfac, zcs_2 ! local scalar 114 CHARACTER (len = 15) :: fieldid 108 REAL(wp) :: za ! local scalar 115 109 REAL(wp), DIMENSION(jpi,jpj) :: closing_net ! net rate at which area is removed (1/s) 116 110 ! ! (ridging ice area - area of new ridges) / dt … … 125 119 IF( nn_timing == 1 ) CALL timing_start('icerdgrft') 126 120 127 IF( kt == nit000 .AND. lwp ) THEN 128 WRITE(numout,*) 129 WRITE(numout,*)'icerdgrft : ice ridging and rafting' 130 WRITE(numout,*)'~~~~~~~~~~' 131 ENDIF 132 !!gm should be: 133 ! IF( kt == nit000 ) THEN 134 ! IF(lwp) WRITE(numout,*) 135 ! IF(lwp) WRITE(numout,*)'icerdgrft : ???' 136 ! IF(lwp) WRITE(numout,*)'~~~~~~~~~~' 137 ! ! 138 ! Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE 139 ! ! 140 !!gm why not adding here zcs_2 computation 141 ! ! 142 ! ENDIF 143 !!gm end 144 121 IF( kt == nit000 ) THEN 122 IF(lwp) WRITE(numout,*) 123 IF(lwp) WRITE(numout,*)'icerdgrft : ice ridging and rafting' 124 IF(lwp) WRITE(numout,*)'~~~~~~~~~~' 125 ! 126 zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE 127 ! 128 ENDIF 145 129 ! ! conservation test 146 130 IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) … … 149 133 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 150 134 !-----------------------------------------------------------------------------! 151 Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE152 zcs_2 = rn_cs * 0.5_wp153 135 ! 154 136 CALL ice_rdgrft_ridgeprep ! prepare ridging … … 175 157 ! (thick, newly ridged ice). 176 158 177 closing_net(ji,jj) = zcs_2* ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )159 closing_net(ji,jj) = rn_cs * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 178 160 179 161 ! 2.2 divu_adv … … 220 202 za = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice 221 203 IF ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN ! would lead to negative ato_i 222 zfac = - ato_i(ji,jj) / za223 204 opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 224 205 ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN ! would lead to ato_i > asum 225 zfac = ( asum(ji,jj) - ato_i(ji,jj) ) / za226 206 opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 227 207 ENDIF … … 238 218 za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 239 219 IF( za > a_i(ji,jj,jl) ) THEN 240 zfac = a_i(ji,jj,jl) / za 241 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 220 closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / za 242 221 ENDIF 243 222 END DO … … 254 233 asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 255 234 256 ! 3.5 Do we keep on iterating ???235 ! 3.5 Do we keep on iterating? 257 236 !-----------------------------------------------------------------------------! 258 237 ! Check whether asum = 1. If not (because the closing and opening … … 316 295 !!---------------------------------------------------------------------! 317 296 INTEGER :: ji, jj, jl ! dummy loop indices 318 REAL(wp) :: Gstari, astari, hrmean, zdummy ! local scalar !!gm DOCTOR norme should start with z !!!!!319 REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: Gsum !Gsum(n) = sum of areas in categories 0 to n297 REAL(wp) :: z1_gstar, z1_astar, zhmean, zdummy ! local scalar 298 REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n 320 299 !------------------------------------------------------------------------------! 321 300 322 Gstari = 1._wp / rn_gstar 323 astari = 1._wp / rn_astar 324 aksum(:,:) = 0._wp 325 athorn(:,:,:) = 0._wp 326 aridge(:,:,:) = 0._wp 327 araft (:,:,:) = 0._wp 301 z1_gstar = 1._wp / rn_gstar 302 z1_astar = 1._wp / rn_astar 328 303 329 304 CALL ice_var_zapsmall ! Zero out categories with very small areas 330 305 331 ! DO jl = 1, jpl ! Ice thickness needed for rafting 332 ! DO jj = 1, jpj 333 ! DO ji = 1, jpi 334 !!gm rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 335 !!gm ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 336 ! IF( a_i(ji,jj,jl) >= epsi20 ) THEN ; ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 337 ! ELSE ; ht_i(ji,jj,jl) = 0._wp 338 ! ENDIF 339 ! END DO 340 ! END DO 341 ! END DO 342 !!gm or even better : 343 ! ! ! Ice thickness needed for rafting 306 ! ! Ice thickness needed for rafting 344 307 WHERE( a_i(:,:,:) >= epsi20 ) ; ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:) 345 308 ELSEWHERE ; ht_i(:,:,:) = 0._wp 346 309 END WHERE 347 !!gm end348 310 349 311 !------------------------------------------------------------------------------! … … 356 318 ! 357 319 ! Compute cumulative thickness distribution function 358 ! Compute the cumulative thickness distribution function Gsum,359 ! where Gsum(n) is the fractional area in categories 0 to n.320 ! Compute the cumulative thickness distribution function zGsum, 321 ! where zGsum(n) is the fractional area in categories 0 to n. 360 322 ! initial value (in h = 0) equals open water area 361 Gsum(:,:,-1) = 0._wp 362 Gsum(:,:,0 ) = ato_i(:,:) 363 ! for each value of h, you have to add ice concentration then 323 zGsum(:,:,-1) = 0._wp 324 zGsum(:,:,0 ) = ato_i(:,:) / asum(:,:) 364 325 DO jl = 1, jpl 365 Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 366 END DO 367 ! 368 ! Normalize the cumulative distribution to 1 369 DO jl = 0, jpl 370 Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:) 326 zGsum(:,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:) ! sum(1:jl) is ok (and not jpl) 371 327 END DO 372 328 … … 388 344 DO jj = 1, jpj 389 345 DO ji = 1, jpi 390 IF ( Gsum(ji,jj,jl) < rn_gstar ) THEN391 athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) -Gsum(ji,jj,jl-1) ) * &392 & ( 2._wp - ( Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari)393 ELSEIF( Gsum(ji,jj,jl-1) < rn_gstar ) THEN394 athorn(ji,jj,jl) = Gstari * ( rn_gstar -Gsum(ji,jj,jl-1) ) * &395 & ( 2._wp - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari)346 IF ( zGsum(ji,jj,jl) < rn_gstar ) THEN 347 athorn(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * & 348 & ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar ) 349 ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN 350 athorn(ji,jj,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jj,jl-1) ) * & 351 & ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar ) * z1_gstar ) 396 352 ELSE 397 353 athorn(ji,jj,jl) = 0._wp … … 403 359 ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007) 404 360 ! 405 zdummy = 1._wp / ( 1._wp - EXP(- astari) ) ! precompute exponential terms usingGsum as a work array361 zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) ) ! precompute exponential terms using zGsum as a work array 406 362 DO jl = -1, jpl 407 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari) * zdummy363 zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy 408 364 END DO 409 365 DO jl = 0, jpl 410 athorn(:,:,jl) = Gsum(:,:,jl-1) -Gsum(:,:,jl)366 athorn(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl) 411 367 END DO 412 368 ! … … 418 374 DO jj = 1, jpj 419 375 DO ji = 1, jpi 420 zdummy = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) 421 aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl) 376 aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * athorn(ji,jj,jl) 422 377 araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl) 423 378 END DO 424 379 END DO 425 380 END DO 426 ELSEIF( ln_ridging .AND. .NOT.ln_rafting ) THEN !- ridging alone 427 DO jl = 1, jpl 428 aridge(:,:,jl) = athorn(:,:,jl) 429 END DO 430 ELSEIF( ln_rafting .AND. .NOT.ln_ridging ) THEN !- rafting alone 431 DO jl = 1, jpl 432 araft(:,:,jl) = athorn(:,:,jl) 433 END DO 381 ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN !- ridging alone 382 aridge(:,:,:) = athorn(:,:,1:jpl) 383 araft (:,:,:) = 0._wp 384 ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone 385 aridge(:,:,:) = 0._wp 386 araft (:,:,:) = athorn(:,:,1:jpl) 387 ELSE !- no ridging & no rafting 388 aridge(:,:,:) = 0._wp 389 araft (:,:,:) = 0._wp 434 390 ENDIF 435 391 … … 441 397 ! 442 398 ! This parameterization is a modified version of Hibler (1980). 443 ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)399 ! The mean ridging thickness, zhmean, is proportional to hi^(0.5) 444 400 ! and for very thick ridging ice must be >= krdgmin*hi 445 401 ! 446 402 ! The minimum ridging thickness, hrmin, is equal to 2*hi 447 403 ! (i.e., rafting) and for very thick ridging ice is 448 ! constrained by hrmin <= ( hrmean + hi)/2.404 ! constrained by hrmin <= (zhmean + hi)/2. 449 405 ! 450 406 ! The maximum ridging thickness, hrmax, is determined by 451 ! hrmean and hrmin.407 ! zhmean and hrmin. 452 408 ! 453 409 ! These modifications have the effect of reducing the ice strength … … 466 422 DO ji = 1, jpi 467 423 IF ( athorn(ji,jj,jl) > 0._wp ) THEN 468 hrmean = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )469 hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) )470 hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl)424 zhmean = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin ) 425 hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) ) 426 hrmax(ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl) 471 427 hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft 472 krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 )428 krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 ) 473 429 ! 474 430 ! Normalization factor : aksum, ensures mass conservation … … 498 454 !!---------------------------------------------------------------------- 499 455 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear 500 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges 501 ! 502 CHARACTER (len=80) :: fieldid ! field identifier !!gm DOCTOR name wrong 456 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area retreats, excluding area of new ridges 503 457 ! 504 458 INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices 505 459 INTEGER :: ij ! horizontal index, combines i and j loops 506 460 INTEGER :: icells ! number of cells with a_i > puny 507 REAL(wp) :: hL, hR, farea ! left and right limits of integration 508 REAL(wp) :: zwfx_snw ! snow mass flux increment 461 REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2 509 462 510 463 INTEGER , DIMENSION(jpij) :: indxi, indxj ! compressed indices 511 REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to n2464 REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2 512 465 513 466 REAL(wp), DIMENSION(jpij) :: afrac ! fraction of category area ridged … … 656 609 ! During the next time step, thermo_rates will determine whether 657 610 ! the ocean cools or new ice grows. 658 zwfx_snw =( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg ) &659 & + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice ! fresh water source for ocean660 661 wfx_snw_dyn(ji,jj) =wfx_snw_dyn(ji,jj) + zwfx_snw611 wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg ) & 612 & + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice ! fresh water source for ocean 613 614 !!! wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + zwfx_snw 662 615 663 616 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg ) & … … 722 675 723 676 ! update jl1 724 e_i 677 e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk) 725 678 726 679 END DO … … 806 759 ! 807 760 INTEGER :: ji,jj, jl ! dummy loop indices 808 INTEGER :: ksmooth ! smoothing the resistance to deformation !!gm not DOCTOR : start with i !!!809 INTEGER :: numts_rm ! number of time steps for the P smoothing !!gm not DOCTOR : start with i !!!761 INTEGER :: ismooth ! smoothing the resistance to deformation 762 INTEGER :: itframe ! number of time steps for the P smoothing 810 763 REAL(wp) :: zp, z1_3 ! local scalars 811 764 REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here … … 813 766 !!---------------------------------------------------------------------- 814 767 815 !------------------------------------------------------------------------------! 816 ! 1) Initialize 817 !------------------------------------------------------------------------------! 818 strength(:,:) = 0._wp 819 820 !------------------------------------------------------------------------------! 821 ! 2) Compute thickness distribution of ridging and ridged ice 822 !------------------------------------------------------------------------------! 823 CALL ice_rdgrft_ridgeprep 824 825 !------------------------------------------------------------------------------! 826 ! 3) Rothrock(1975)'s method 827 !------------------------------------------------------------------------------! 828 IF( kstrngth == 1 ) THEN 768 ! !--------------------------------------------------! 769 CALL ice_rdgrft_ridgeprep ! Thickness distribution of ridging and ridged ice ! 770 ! !--------------------------------------------------! 771 772 ! !--------------------------------------------------! 773 IF( kstrngth == 1 ) THEN ! Ice strength => Rothrock (1975) method ! 774 ! !--------------------------------------------------! 829 775 z1_3 = 1._wp / 3._wp 830 776 DO jl = 1, jpl 831 DO jj= 1, jpj 832 DO ji = 1, jpi 833 ! 834 IF( athorn(ji,jj,jl) > 0._wp ) THEN 835 !---------------------------- 836 ! PE loss from deforming ice 837 !---------------------------- 838 strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl) 839 840 !-------------------------- 841 ! PE gain from rafting ice 842 !-------------------------- 843 strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl) 844 845 !---------------------------- 846 ! PE gain from ridging ice 847 !---------------------------- 848 strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 * & 849 & ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) + & 850 & hrmin(ji,jj,jl) * hrmin(ji,jj,jl) + & 851 & hrmax(ji,jj,jl) * hrmin(ji,jj,jl) ) 852 !!(a**3-b**3)/(a-b) = a*a+ab+b*b 853 ENDIF 854 ! 855 END DO 856 END DO 857 END DO 858 859 strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1) 860 ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 861 ksmooth = 1 862 863 !------------------------------------------------------------------------------! 864 ! 4) Hibler (1979)' method 865 !------------------------------------------------------------------------------! 866 ELSE ! kstrngth ne 1: Hibler (1979) form 867 ! 777 WHERE( athorn(:,:,jl) > 0._wp ) 778 strength(:,:) = - athorn(:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl) & ! PE loss from deforming ice 779 & + 2._wp * araft (:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl) & ! PE gain from rafting ice 780 & + aridge(:,:,jl) * krdg(:,:,jl) * z1_3 * & ! PE gain from ridging ice 781 & ( hrmax(:,:,jl) * hrmax(:,:,jl) + & 782 & hrmin(:,:,jl) * hrmin(:,:,jl) + & 783 & hrmax(:,:,jl) * hrmin(:,:,jl) ) 784 ELSEWHERE 785 strength(:,:) = 0._wp 786 END WHERE 787 END DO 788 strength(:,:) = rn_pe_rdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1) 789 ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 790 ismooth = 1 791 ! !--------------------------------------------------! 792 ELSE ! Ice strength => Hibler (1979) method ! 793 ! !--------------------------------------------------! 868 794 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) * tmask(:,:,1) 869 795 ! 870 ksmooth = 1 871 ! 872 ENDIF ! kstrngth 873 ! 874 !------------------------------------------------------------------------------! 875 ! 5) Impact of brine volume 876 !------------------------------------------------------------------------------! 796 ismooth = 1 797 ! 798 ENDIF 799 ! !--------------------------------------------------! 800 IF( ln_icestr_bvf ) THEN ! Impact of brine volume ! 801 ! !--------------------------------------------------! 877 802 ! CAN BE REMOVED 878 IF( ln_icestr_bvf ) THEN879 803 DO jj = 1, jpj 880 804 DO ji = 1, jpi … … 883 807 END DO 884 808 ENDIF 885 ! 886 !------------------------------------------------------------------------------! 887 ! 6) Smoothing ice strength 888 !------------------------------------------------------------------------------! 889 SELECT CASE( ksmooth ) 890 ! !------------------- 891 CASE( 1 ) ! Spatial smoothing 892 ! !------------------- 809 ! !--------------------------------------------------! 810 SELECT CASE( ismooth ) ! Smoothing ice strength ! 811 ! !--------------------------------------------------! 812 ! 813 CASE( 1 ) !--- Spatial smoothing 893 814 DO jj = 2, jpjm1 894 815 DO ji = 2, jpim1 … … 903 824 END DO 904 825 END DO 905 826 906 827 DO jj = 2, jpjm1 907 828 DO ji = 2, jpim1 … … 911 832 CALL lbc_lnk( strength, 'T', 1. ) 912 833 ! 913 ! !-------------------- 914 CASE( 2 ) ! Temporal smoothing 915 ! !-------------------- 834 CASE( 2 ) !--- Temporal smoothing 916 835 IF ( kt_ice == nit000 ) THEN 917 836 zstrp1(:,:) = 0._wp … … 922 841 DO ji = 2, jpim1 923 842 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN 924 numts_rm= 1 ! number of time steps for the running mean925 IF ( zstrp1(ji,jj) > 0._wp ) numts_rm = numts_rm+ 1926 IF ( zstrp2(ji,jj) > 0._wp ) numts_rm = numts_rm+ 1927 zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / numts_rm843 itframe = 1 ! number of time steps for the running mean 844 IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1 845 IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1 846 zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe 928 847 zstrp2 (ji,jj) = zstrp1 (ji,jj) 929 848 zstrp1 (ji,jj) = strength(ji,jj) … … 932 851 END DO 933 852 END DO 934 CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions853 CALL lbc_lnk( strength, 'T', 1. ) 935 854 ! 936 855 END SELECT … … 970 889 IF (lwp) THEN ! control print 971 890 WRITE(numout,*) 972 WRITE(numout,*)'ice_rdgrft_init : ice parameters for mechanical ice redistribution ' 973 WRITE(numout,*)'~~~~~~~~~~~~~~~' 974 WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs 975 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun 976 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar 977 WRITE(numout,*)' Equivalent to G* for an exponential part function rn_astar = ', rn_astar 978 WRITE(numout,*)' Ridging of ice sheets or not ln_ridging = ', ln_ridging 979 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar 980 WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg 981 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg 982 WRITE(numout,*)' Fraction of pond volume conserved during ridging rn_fpondrdg = ', rn_fpondrdg 983 WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting 984 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft 985 WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft 986 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft 987 WRITE(numout,*)' Fraction of pond volume conserved during rafting rn_fpondrft = ', rn_fpondrft 891 WRITE(numout,*) 'ice_rdgrft_init : ice parameters for mechanical ice redistribution ' 892 WRITE(numout,*) '~~~~~~~~~~~~~~~' 893 WRITE(numout,*) ' Namelist namiceitdme' 894 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs 895 WRITE(numout,*) ' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun 896 WRITE(numout,*) ' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar 897 WRITE(numout,*) ' Equivalent to G* for an exponential part function rn_astar = ', rn_astar 898 WRITE(numout,*) ' Ridging of ice sheets or not ln_ridging = ', ln_ridging 899 WRITE(numout,*) ' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar 900 WRITE(numout,*) ' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg 901 WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg 902 WRITE(numout,*) ' Fraction of pond volume conserved during ridging rn_fpondrdg = ', rn_fpondrdg 903 WRITE(numout,*) ' Rafting of ice sheets or not ln_rafting = ', ln_rafting 904 WRITE(numout,*) ' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft 905 WRITE(numout,*) ' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft 906 WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft 907 WRITE(numout,*) ' Fraction of pond volume conserved during rafting rn_fpondrft = ', rn_fpondrft 988 908 ENDIF 989 909 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
r8486 r8498 31 31 USE prtctl ! Print control 32 32 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 33 USE iom ! 33 34 #if defined key_agrif 34 35 USE agrif_lim3_interp … … 111 112 INTEGER :: ji, jj ! dummy loop indices 112 113 INTEGER :: jter ! local integers 113 CHARACTER (len=50) :: charout114 114 115 115 REAL(wp) :: zrhoco ! rau0 * rn_cio … … 121 121 REAL(wp) :: zTauO, zTauB, zTauE, zvel ! temporary scalars 122 122 123 REAL(wp) :: zsig1, zsig2 ! internal ice stress124 123 REAL(wp) :: zresm ! Maximal error on ice velocity 125 124 REAL(wp) :: zintb, zintn ! dummy argument 126 125 REAL(wp) :: zfac_x, zfac_y 126 REAL(wp) :: zshear, zdum1, zdum2 127 127 128 128 REAL(wp), DIMENSION(jpi,jpj) :: z1_e1t0, z1_e2t0 ! scale factors … … 152 152 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 153 153 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity 154 !! --- diags 155 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zswi, zsig1, zsig2, zsig3 156 !! --- SIMIP diags 157 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig1 ! Average normal stress in sea ice 158 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig2 ! Maximum shear stress in sea ice 159 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_dssh_dx ! X-direction sea-surface tilt term (N/m2) 160 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_dssh_dy ! X-direction sea-surface tilt term (N/m2) 161 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_corstrx ! X-direction coriolis stress (N/m2) 162 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_corstry ! Y-direction coriolis stress (N/m2) 163 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_intstrx ! X-direction internal stress (N/m2) 164 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_intstry ! Y-direction internal stress (N/m2) 165 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_utau_oi ! X-direction ocean-ice stress 166 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_vtau_oi ! Y-direction ocean-ice stress 167 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 168 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) 169 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s) 170 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) 171 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s) 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s) 154 173 !!------------------------------------------------------------------- 155 174 … … 671 690 672 691 !------------------------------------------------------------------------------! 673 ! 5) SIMIP diagnostics 674 !------------------------------------------------------------------------------! 675 676 !!gm encapsulate with a flag (iom_use of the variable or better a flag defined one for all testing if one of the 677 !! diag is output. NB the diag_... are should only be ALLOCATED if the flag is true ! 678 679 DO jj = 2, jpjm1 680 DO ji = 2, jpim1 681 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 682 683 ! Stress tensor invariants (normal and shear stress N/m) 684 diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch ! normal stress 685 diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch ! shear stress 686 687 ! Stress terms of the momentum equation (N/m2) 688 diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch ! sea surface slope stress term 689 diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 690 691 diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch ! Coriolis stress term 692 diag_corstry(ji,jj) = zCory(ji,jj) * rswitch 693 694 diag_intstrx(ji,jj) = zfU(ji,jj) * rswitch ! internal stress term 695 diag_intstry(ji,jj) = zfV(ji,jj) * rswitch 696 697 diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch ! oceanic stress 698 diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 699 700 ! 2D ice mass, snow mass, area transport arrays (X, Y) 701 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 702 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 703 704 diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 705 diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' 706 707 diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 708 diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' 709 710 diag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component 711 diag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 712 713 END DO 714 END DO 715 716 CALL lbc_lnk_multi( diag_sig1 , 'T', 1., diag_sig2 , 'T', 1., & 717 & diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1., & 718 & diag_corstrx, 'U', -1., diag_corstry, 'V', -1., & 719 & diag_intstrx, 'U', -1., diag_intstry, 'V', -1. ) 720 721 CALL lbc_lnk_multi( diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1. ) 722 723 CALL lbc_lnk_multi( diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1., & 724 & diag_xatrp , 'U', -1., diag_ymtrp_ice, 'V', -1., & 725 & diag_ymtrp_snw, 'V', -1., diag_yatrp , 'V', -1. ) 726 727 ! 728 !------------------------------------------------------------------------------! 729 ! 6) Control prints of residual and charge ellipse 730 !------------------------------------------------------------------------------! 731 ! 732 ! print the residual for convergence 733 IF(ln_ctl) THEN 734 WRITE(charout,FMT="('ice_rhg_evp : res =',D23.16, ' iter =',I4)") zresm, jter 735 CALL prt_ctl_info(charout) 736 CALL prt_ctl(tab2d_1=u_ice, clinfo1=' ice_rhg_evp : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 692 ! 5) diagnostics 693 !------------------------------------------------------------------------------! 694 ! --- ellipse --- ! 695 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN 737 696 ! 738 ! print charge ellipse 739 ! This can be desactivated once the user is sure that the stress state 740 ! lie on the charge ellipse. See Bouillon et al. (2008) for more details 741 IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN 742 WRITE(charout,FMT="('ice_rhg_evp :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. ' 743 CALL prt_ctl_info(charout) 744 DO jj = 2, jpjm1 745 DO ji = 2, jpim1 746 IF (strength(ji,jj) > 1.0) THEN 747 zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 748 zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 749 WRITE(charout,FMT="('ice_rhg_evp :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 750 CALL prt_ctl_info(charout) 751 ENDIF 752 END DO 753 END DO 754 WRITE(charout,FMT="('ice_rhg_evp :', I4, I6, I1, I1, A10)") 2000, kt_ice, 0, 0, ' ch. ell. ' 755 CALL prt_ctl_info(charout) 756 ENDIF 697 ALLOCATE( zswi(jpi,jpj) , zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 698 ! 699 DO jj = 1, jpj 700 DO ji = 1, jpi 701 zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 702 END DO 703 END DO 704 705 DO jj = 2, jpjm1 706 DO ji = 2, jpim1 707 zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji ,jj-1) * stress12_i(ji ,jj-1) + & ! stress12_i at T-point 708 & zswi(ji ,jj) * stress12_i(ji ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) ) & 709 & / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 710 711 zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 712 713 zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 714 715 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 716 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 717 !! zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 718 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 719 zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 720 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 721 zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 722 END DO 723 END DO 724 CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 725 ! 726 IF( iom_use('isig1') ) CALL iom_put( "isig1" , zsig1 ) 727 IF( iom_use('isig2') ) CALL iom_put( "isig2" , zsig2 ) 728 IF( iom_use('isig3') ) CALL iom_put( "isig3" , zsig3 ) 729 ! 730 DEALLOCATE( zswi , zsig1 , zsig2 , zsig3 ) 731 ENDIF 732 733 ! --- SIMIP --- ! 734 IF ( iom_use( 'normstr' ) .OR. iom_use( 'sheastr' ) .OR. iom_use( 'dssh_dx' ) .OR. iom_use( 'dssh_dy' ) .OR. & 735 & iom_use( 'corstrx' ) .OR. iom_use( 'corstry' ) .OR. iom_use( 'intstrx' ) .OR. iom_use( 'intstry' ) .OR. & 736 & iom_use( 'utau_oi' ) .OR. iom_use( 'vtau_oi' ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. & 737 & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp' ) .OR. iom_use( 'yatrp' ) ) THEN 738 739 ALLOCATE( zdiag_sig1 (jpi,jpj) , zdiag_sig2 (jpi,jpj) , zdiag_dssh_dx (jpi,jpj) , zdiag_dssh_dy (jpi,jpj) , & 740 & zdiag_corstrx (jpi,jpj) , zdiag_corstry (jpi,jpj) , zdiag_intstrx (jpi,jpj) , zdiag_intstry (jpi,jpj) , & 741 & zdiag_utau_oi (jpi,jpj) , zdiag_vtau_oi (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & 742 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp (jpi,jpj) , zdiag_yatrp (jpi,jpj) ) 743 744 DO jj = 2, jpjm1 745 DO ji = 2, jpim1 746 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 747 748 ! Stress tensor invariants (normal and shear stress N/m) 749 zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch ! normal stress 750 zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch ! shear stress 751 752 ! Stress terms of the momentum equation (N/m2) 753 zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch ! sea surface slope stress term 754 zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 755 756 zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch ! Coriolis stress term 757 zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch 758 759 zdiag_intstrx(ji,jj) = zfU(ji,jj) * rswitch ! internal stress term 760 zdiag_intstry(ji,jj) = zfV(ji,jj) * rswitch 761 762 zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch ! oceanic stress 763 zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 764 765 ! 2D ice mass, snow mass, area transport arrays (X, Y) 766 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 767 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 768 769 zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 770 zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' 771 772 zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 773 zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' 774 775 zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component 776 zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 777 778 END DO 779 END DO 780 781 CALL lbc_lnk_multi( zdiag_sig1 , 'T', 1., zdiag_sig2 , 'T', 1., & 782 & zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1., & 783 & zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1., & 784 & zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1. ) 785 786 CALL lbc_lnk_multi( zdiag_utau_oi, 'U', -1., zdiag_vtau_oi, 'V', -1. ) 787 788 CALL lbc_lnk_multi( zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1., & 789 & zdiag_xatrp , 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 790 & zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp , 'V', -1. ) 791 792 IF ( iom_use( 'normstr' ) ) CALL iom_put( 'normstr' , zdiag_sig1(:,:) ) ! Normal stress 793 IF ( iom_use( 'sheastr' ) ) CALL iom_put( 'sheastr' , zdiag_sig2(:,:) ) ! Shear stress 794 IF ( iom_use( 'dssh_dx' ) ) CALL iom_put( 'dssh_dx' , zdiag_dssh_dx(:,:) ) ! Sea-surface tilt term in force balance (x) 795 IF ( iom_use( 'dssh_dy' ) ) CALL iom_put( 'dssh_dy' , zdiag_dssh_dy(:,:) ) ! Sea-surface tilt term in force balance (y) 796 IF ( iom_use( 'corstrx' ) ) CALL iom_put( 'corstrx' , zdiag_corstrx(:,:) ) ! Coriolis force term in force balance (x) 797 IF ( iom_use( 'corstry' ) ) CALL iom_put( 'corstry' , zdiag_corstry(:,:) ) ! Coriolis force term in force balance (y) 798 IF ( iom_use( 'intstrx' ) ) CALL iom_put( 'intstrx' , zdiag_intstrx(:,:) ) ! Internal force term in force balance (x) 799 IF ( iom_use( 'intstry' ) ) CALL iom_put( 'intstry' , zdiag_intstry(:,:) ) ! Internal force term in force balance (y) 800 IF ( iom_use( 'utau_oi' ) ) CALL iom_put( 'utau_oi' , zdiag_utau_oi(:,:) ) ! Ocean stress term in force balance (x) 801 IF ( iom_use( 'vtau_oi' ) ) CALL iom_put( 'vtau_oi' , zdiag_vtau_oi(:,:) ) ! Ocean stress term in force balance (y) 802 IF ( iom_use( 'xmtrpice' ) ) CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice(:,:) ) ! X-component of sea-ice mass transport (kg/s) 803 IF ( iom_use( 'ymtrpice' ) ) CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice(:,:) ) ! Y-component of sea-ice mass transport 804 IF ( iom_use( 'xmtrpsnw' ) ) CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw(:,:) ) ! X-component of snow mass transport (kg/s) 805 IF ( iom_use( 'ymtrpsnw' ) ) CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw(:,:) ) ! Y-component of snow mass transport 806 IF ( iom_use( 'xatrp' ) ) CALL iom_put( 'xatrp' , zdiag_xatrp(:,:) ) ! X-component of ice area transport 807 IF ( iom_use( 'yatrp' ) ) CALL iom_put( 'yatrp' , zdiag_yatrp(:,:) ) ! Y-component of ice area transport 808 809 DEALLOCATE( zdiag_sig1 , zdiag_sig2 , zdiag_dssh_dx , zdiag_dssh_dy , & 810 & zdiag_corstrx , zdiag_corstry , zdiag_intstrx , zdiag_intstry , & 811 & zdiag_utau_oi , zdiag_vtau_oi , zdiag_xmtrp_ice , zdiag_ymtrp_ice , & 812 & zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) 813 757 814 ENDIF 758 815 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
r8486 r8498 122 122 ! !-----------------------! 123 123 ! 124 kt_ice = kt !Ice model time step124 kt_ice = kt ! -- Ice model time step 125 125 ! 126 126 # if defined key_agrif 127 127 IF( .NOT. Agrif_Root() ) lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1 128 128 # endif 129 130 ! ! mean surface ocean current at ice velocity point 131 u_oce(:,:) = ssu_m(:,:) * umask(:,:,1) 132 v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1) 133 !!gm the umask, vmask above is useless as ssu_m, ssv_m are build from masked un,vn (used t be here for B-grid sea-ice) 134 135 ! ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 136 CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) 129 u_oce(:,:) = ssu_m(:,:) ! -- mean surface ocean current 130 v_oce(:,:) = ssv_m(:,:) 131 ! 132 CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) ! -- freezing temperature [Kelvin] (set to rt0 over land) 137 133 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 138 134 ! 139 CALL ice_bef ! Store previousice values135 CALL store_fields ! -- Store now ice values 140 136 ! 141 137 !------------------------------------------------! … … 179 175 CALL ice_var_glo2eqv ! ht_i and ht_s for ice albedo calculation 180 176 CALL ice_var_agg(1) ! at_i for coupling 181 CALL ice_bef ! Store previousice values177 CALL store_fields ! Store now ice values 182 178 183 179 !------------------------------------------------------! … … 223 219 !! IF( .NOT. Agrif_Root() ) CALL Agrif_ParentGrid_To_ChildGrid() 224 220 !!# endif 225 IF( ln_limdiahsb ) 226 ! 227 221 IF( ln_limdiahsb ) CALL ice_dia( kt ) ! -- Diagnostics and outputs 222 ! 223 CALL ice_wri( 1 ) ! -- Ice outputs 228 224 ! 229 225 IF( kt == nit000 .AND. ln_rstart ) & !!gm I don't understand the ln_rstart, if needed, add a comment, please 230 & 231 ! 232 IF( lrst_ice ) 233 ! 234 IF( ln_limctl ) 226 & CALL iom_close( numrir ) ! close input ice restart file 227 ! 228 IF( lrst_ice ) CALL ice_rst_write( kt ) ! -- Ice restart file 229 ! 230 IF( ln_limctl ) CALL ice_ctl( kt ) ! alerts in case of model crash 235 231 ! 236 232 ENDIF ! End sea-ice time step only … … 241 237 ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere 242 238 ! using before instantaneous surf. currents 243 IF( ln_limdyn ) 239 IF( ln_limdyn ) CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) ) 244 240 !!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!! 245 241 ! … … 476 472 477 473 478 SUBROUTINE ice_bef479 !!---------------------------------------------------------------------- 480 !! *** ROUTINE ice_bef***474 SUBROUTINE store_fields 475 !!---------------------------------------------------------------------- 476 !! *** ROUTINE store_fields *** 481 477 !! 482 478 !! ** purpose : store ice variables at "before" time step … … 525 521 CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. ) 526 522 ! 527 END SUBROUTINE ice_bef523 END SUBROUTINE store_fields 528 524 529 525 … … 570 566 ! 571 567 afx_tot(ji,jj) = 0._wp ; 572 afx_dyn(ji,jj) = 0._wp ; afx_thd(ji,jj) = 0._wp573 568 ! 574 569 diag_heat(ji,jj) = 0._wp ; diag_smvi(ji,jj) = 0._wp -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
r8486 r8498 193 193 ! Second step in icethd_dh : heat remaining if total melt (zq_rema) 194 194 ! Third step in iceupdate.F90 : heat from ice-ocean mass exchange (zf_mass) + solar 195 DO jj = 1, jpj 196 DO ji = 1, jpi 197 hfx_out(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) & ! Non solar heat flux received by the ocean 198 & - qlead(ji,jj) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation 199 & - at_i(ji,jj) * fhtur(ji,jj) & ! heat flux taken by turbulence 200 & - at_i(ji,jj) * fhld(ji,jj) ! heat flux taken during bottom growth/melt 201 ! (fhld should be 0 while bott growth) 202 END DO 203 END DO 204 195 hfx_out(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) & ! Non solar heat flux received by the ocean 196 & - qlead(:,:) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation 197 & - at_i (:,:) * fhtur(:,:) & ! heat flux taken by turbulence 198 & - at_i (:,:) * fhld(:,:) ! heat flux taken during bottom growth/melt 199 ! (fhld should be 0 while bott growth) 205 200 !-------------------------------------------------------------------------------------------! 206 201 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories … … 288 283 !!------------------------------------------------------------------- 289 284 INTEGER :: ji, jk ! dummy loop indices 290 REAL(wp) :: ztmelts, z aaa, zbbb, zccc, zdiscrim! local scalar285 REAL(wp) :: ztmelts, z1_2cp, zbbb, zccc ! local scalar 291 286 !!------------------------------------------------------------------- 292 287 ! Recover ice temperature 288 z1_2cp = 1._wp / ( 2._wp * cpic ) 293 289 DO jk = 1, nlay_i 294 290 DO ji = 1, nidx 295 ztmelts = -tmut * s_i_1d(ji,jk) + rt0291 ztmelts = -tmut * s_i_1d(ji,jk) 296 292 ! Conversion q(S,T) -> T (second order equation) 297 zaaa = cpic 298 zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + e_i_1d(ji,jk) * r1_rhoic - lfus 299 zccc = lfus * ( ztmelts - rt0 ) 300 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 301 t_i_1d(ji,jk) = rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 293 zbbb = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus 294 zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) ) 295 t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * z1_2cp 302 296 303 297 ! mask temperature 304 rswitch = 305 t_i_1d(ji,jk) = 298 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 299 t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0 306 300 END DO 307 301 END DO -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90
r8486 r8498 32 32 CONTAINS 33 33 34 !!gm even comment line of more than 130 character may cause compilation errors35 !!gm ===>>> reformat the text36 34 SUBROUTINE ice_thd_da 37 35 !!------------------------------------------------------------------- … … 46 44 !! W = m1 * (Tw -Tf)**m2 --- originally from Josberger 1979 --- 47 45 !! (Tw - Tf) = elevation of water temp above freezing 48 !! m1 and m2 = (1.6e-6 , 1.36) best fit from field experiment near the coast of Prince Patrick Island (Perovich 1983) => static ice 49 !! m1 and m2 = (3.0e-6 , 1.36) best fit from MIZEX 84 experiment (Maykut and Perovich 1987) => moving ice 46 !! m1 and m2 = (1.6e-6 , 1.36) best fit from field experiment near the coast of Prince Patrick Island 47 !! (Perovich 1983) => static ice 48 !! m1 and m2 = (3.0e-6 , 1.36) best fit from MIZEX 84 experiment 49 !! (Maykut and Perovich 1987) => moving ice 50 50 !! 51 51 !! P = N * pi * D --- from Rothrock and Thorndike 1984 --- … … 59 59 !! Astar = 1 / ( 1 - (Dmin/Dmax)**(1/beta) ) 60 60 !! Dmin = minimum floe diameter (recommended to be 8m +- 20%) 61 !! Dmax = maximum floe diameter (recommended to be 300m, but it does not impact melting much except for Dmax<100m) 61 !! Dmax = maximum floe diameter (recommended to be 300m, 62 !! but it does not impact melting much except for Dmax<100m) 62 63 !! beta = 1.0 +-20% (recommended value) 63 64 !! = 0.3 best fit for western Fram Strait and Antarctica … … 113 114 ! --- Calculate reduction of total sea ice concentration --- ! 114 115 zdfloe = rn_dmin * ( zastar / ( zastar - at_i_1d(ji) ) )**rn_beta ! Mean floe caliper diameter [m] 115 zperi = at_i_1d(ji) * rpi / ( zcs * zdfloe ) ! Mean perimeter of the floe = N*pi*D = (A/cs*D^2)*pi*D [m.m-2] 116 ! 117 zperi = at_i_1d(ji) * rpi / ( zcs * zdfloe ) ! Mean perimeter of the floe [m.m-2] 118 ! ! = N*pi*D = (A/cs*D^2)*pi*D 116 119 zwlat = zm1 * ( MAX( 0._wp, sst_1d(ji) - ( t_bo_1d(ji) - rt0 ) ) )**zm2 ! Melt speed rate [m/s] 117 120 ! 118 121 zda_tot(ji) = MIN( zwlat * zperi * rdt_ice, at_i_1d(ji) ) ! sea ice concentration decrease (>0) 119 122 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
r8486 r8498 102 102 REAL(wp), DIMENSION(jpij) :: zsnw ! distribution of snow after wind blowing 103 103 104 REAL(wp) :: zswitch_sal104 REAL(wp) :: zswitch_sal 105 105 106 106 INTEGER :: num_iter_max ! Heat conservation -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90
r8486 r8498 53 53 !!--------------------------------------------------------------------- 54 54 55 !--------------------------------------------------------------------| 56 ! 1) salinity constant in time | 57 !--------------------------------------------------------------------| 58 ! do nothing 59 60 !----------------------------------------------------------------------| 61 ! 2) salinity varying in time | 62 !----------------------------------------------------------------------| 63 IF( nn_icesal == 2 ) THEN 64 55 SELECT CASE ( nn_icesal ) 56 ! 57 ! !---------------------------------------------! 58 CASE( 2 ) ! time varying salinity with linear profile ! 59 ! !---------------------------------------------! 65 60 DO ji = 1, nidx 66 61 … … 96 91 CALL ice_var_salprof1d 97 92 ! 98 ENDIF 99 100 !------------------------------------------------------------------------------| 101 ! 3) vertical profile of salinity, constant in time | 102 !------------------------------------------------------------------------------| 103 IF( nn_icesal == 3 ) CALL ice_var_salprof1d 93 ! !---------------------------------------------! 94 CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile(mean = 2.30) 95 ! !---------------------------------------------! 96 CALL ice_var_salprof1d 104 97 ! 98 END SELECT 99 ! 105 100 END SUBROUTINE ice_thd_sal 106 101 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
r8486 r8498 86 86 !! 87 87 !! ** Purpose : Update the surface ocean boundary condition for heat 88 !! salt and mass over areas where sea-ice is non-zero88 !! salt and mass over areas where sea-ice is non-zero 89 89 !! 90 90 !! ** Action : - computes the heat and freshwater/salt fluxes 91 !! at the ice-ocean interface.91 !! at the ice-ocean interface. 92 92 !! - Update the ocean sbc 93 93 !! … … 133 133 DO ji = 1, jpi 134 134 135 !------------------------------------------!136 ! heat flux at the ocean surface !137 !------------------------------------------!138 135 ! Solar heat flux reaching the ocean = zqsr (W.m-2) 139 136 !--------------------------------------------------- 140 zqsr = qsr_tot(ji,jj) 141 DO jl = 1, jpl 142 zqsr = zqsr - a_i_b(ji,jj,jl) * ( qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) ) 143 END DO 144 !!gm why not like almost everywhere else : 145 !!gm zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) ) 137 zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) ) ) 146 138 147 139 ! Total heat flux reaching the ocean = hfx_out (W.m-2) … … 204 196 END DO 205 197 206 !-----------------------------------------------! 207 ! Storing the transmitted variables ! 208 !-----------------------------------------------! 198 ! Storing the transmitted variables 199 !---------------------------------- 209 200 fr_i (:,:) = at_i(:,:) ! Sea-ice fraction 210 201 tn_ice(:,:,:) = t_su(:,:,:) ! Ice surface temperature 211 202 212 !------------------------------------------------------------------------! 213 ! Snow/ice albedo (only if sent to coupler, useless in forced mode) ! 214 !------------------------------------------------------------------------! 203 ! Snow/ice albedo (only if sent to coupler, useless in forced mode) 204 !------------------------------------------------------------------ 215 205 CALL ice_alb( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 216 206 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
r8496 r8498 157 157 !!------------------------------------------------------------------ 158 158 INTEGER :: ji, jj, jk, jl ! dummy loop indices 159 REAL(wp) :: ze_i, z1_ CpR, zdiscrim, zaaa, z1_2aaa! local scalars160 REAL(wp) :: ze_s, z L_Cp , ztmelts , zbbb, zccc! - -161 REAL(wp) :: zhmax, z1_zhmax, zsm_i , zcpMcpic, zt_i! - -162 REAL(wp) :: zlay_i, zlay_s ! - -159 REAL(wp) :: ze_i, z1_cp, z1_2cp ! local scalars 160 REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - - 161 REAL(wp) :: zhmax, z1_zhmax, zsm_i ! - - 162 REAL(wp) :: zlay_i, zlay_s ! - - 163 163 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i 164 164 !!------------------------------------------------------------------ … … 181 181 ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) !--- ice thickness 182 182 183 183 zhmax = hi_max(jpl) 184 184 z1_zhmax = 1._wp / hi_max(jpl) 185 185 WHERE( ht_i(:,:,jpl) > zhmax ) !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area … … 205 205 !------------------- 206 206 zlay_i = REAL( nlay_i , wp ) ! number of layers 207 zaaa = cpic ! Conversion q(S,T) -> T (second order equation) 208 z1_2aaa = 1._wp / ( 2._wp * zaaa ) 209 zcpMcpic = rcp - cpic 207 z1_2cp = 1._wp / ( 2._wp * cpic ) 210 208 DO jl = 1, jpl 211 209 DO jk = 1, nlay_i … … 214 212 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 215 213 ! 216 ze_i = e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 217 ztmelts = - s_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C] 218 ! 219 zbbb = zcpMcpic * ztmelts + ze_i * r1_rhoic - lfus 220 zccc = lfus * ztmelts 221 zdiscrim = SQRT( MAX( zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 222 zt_i = - ( zbbb + zdiscrim ) * z1_2aaa ! [C] 223 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( zt_i , ztmelts ) ) + rt0 ! [K] with bounds: -100 < zt_i < ztmelts 214 ze_i = e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 215 ztmelts = - s_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C] 216 ! Conversion q(S,T) -> T (second order equation) 217 zbbb = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus 218 zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) ) 219 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * z1_2cp , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts 224 220 ! 225 221 ELSE !--- no ice 226 t_i(ji,jj,jk,jl) = 222 t_i(ji,jj,jk,jl) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C) 227 223 !!clem: I think we should impose rt0 instead 228 224 ENDIF … … 235 231 ! Snow temperature [K] (with a minimum value (rt0 - 100.) imposed everywhere) 236 232 !-------------------- 237 z1_CpR = 1._wp / ( cpic * rhosn )238 zL_Cp = lfus / cpic239 233 zlay_s = REAL( nlay_s , wp ) 234 z1_cp = 1._wp / cpic 240 235 DO jk = 1, nlay_s 241 236 WHERE( v_s(:,:,:) > epsi20 ) !--- icy area 242 t_s(:,:,jk,:) = MAX( -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp ) ) + rt0237 t_s(:,:,jk,:) = MAX( -100._wp , MIN( z1_cp * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0 243 238 ELSEWHERE !--- no ice 244 239 t_s(:,:,jk,:) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C) 245 240 END WHERE 246 241 END DO 247 !!gm perhaps better like this (?) : (if this compile, yes! one test instead of nlay_s tests)248 ! WHERE( v_s(:,:,:) > epsi20 ) !--- icy area249 ! DO jk = 1, nlay_s250 ! t_s(:,:,jk,:) = MAX( -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp ) ) + rt0251 ! END DO252 ! ELSEWHERE !--- no ice253 ! DO jk = 1, nlay_s254 ! t_s(:,:,jk,:) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C)255 ! END DO256 ! END WHERE257 !!gm end better ?258 242 259 243 ! integrated values -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90
r8488 r8498 54 54 INTEGER :: ji, jj, jk, jl ! dummy loop indices 55 55 REAL(wp) :: z2da, z2db, ztmp, zrho1, zrho2, zmiss_val 56 REAL(wp) :: zs12, zshear57 56 REAL(wp), DIMENSION(jpi,jpj) :: z2d, zswi, zmiss ! 2D workspace 58 57 REAL(wp), DIMENSION(jpi,jpj) :: zfb ! ice freeboard 59 58 REAL(wp), DIMENSION(jpi,jpj) :: zamask, zamask15 ! 15% concentration mask 60 REAL(wp), DIMENSION(jpi,jpj) :: zsig1, zsig2, zsig361 59 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zswi2, zmiss2 62 60 ! … … 196 194 CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp ) ! "excess" of sublimation sent to ocean 197 195 198 CALL iom_put( "afxtot" , afx_tot ) ! concentration tendency (total)199 CALL iom_put( "afxdyn" , afx_dyn ) ! concentration tendency (dynamics)200 CALL iom_put( "afxthd" , afx_thd ) ! concentration tendency (thermo)201 202 196 CALL iom_put ('hfxthd' , hfx_thd(:,:) ) ! 203 197 CALL iom_put ('hfxdyn' , hfx_dyn(:,:) ) ! … … 219 213 CALL iom_put ('hfxspr' , hfx_spr(:,:) ) ! Heat content of snow precip 220 214 221 !!gm ====>>>>> THIS should be moved in icerhg_evp (generalize this everywhere it is possible and logic...)222 ! specific outputs for EVP rheology223 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN224 zsig1(:,:) = 0._wp; zsig2(:,:) = 0._wp; zsig3(:,:) = 0._wp;225 DO jj = 2, jpjm1226 DO ji = 2, jpim1227 zs12 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji ,jj-1) * stress12_i(ji ,jj-1) + & ! stress12_i at T-point228 & zswi(ji ,jj) * stress12_i(ji ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) ) &229 & / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )230 231 zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zs12 * zs12 ) ! shear stress232 233 z2da = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )234 235 !! zsig1(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)236 !! zsig2(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)237 !! zsig3(ji,jj) = z2da**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress238 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))239 zsig1(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015240 zsig2(ji,jj) = 0.5_wp * z2da * ( zshear ) ! shear stress241 zsig3(ji,jj) = z2da**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )242 END DO243 END DO244 CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )245 CALL iom_put( "isig1" , zsig1 )246 CALL iom_put( "isig2" , zsig2 )247 CALL iom_put( "isig3" , zsig3 )248 ENDIF249 !!gm <<<<<<======= end250 251 215 !!gm ====>>>>> THIS should be moved where at_ip, vt_ip are computed fro the last time in the time-step (limmpd) 252 216 ! MV MP 2016 … … 339 303 IF ( iom_use( "vice_mv" ) ) CALL iom_put( "vice_mv" , v_ice(:,:) * zswi(:,:) + zmiss(:,:) ) ! ice velocity v component 340 304 341 IF ( iom_use( "xmtrpice" ) ) CALL iom_put( "xmtrpice" , diag_xmtrp_ice(:,:) ) ! X-component of sea-ice mass transport (kg/s)342 IF ( iom_use( "ymtrpice" ) ) CALL iom_put( "ymtrpice" , diag_ymtrp_ice(:,:) ) ! Y-component of sea-ice mass transport343 344 IF ( iom_use( "xmtrpsnw" ) ) CALL iom_put( "xmtrpsnw" , diag_xmtrp_snw(:,:) ) ! X-component of snow mass transport (kg/s)345 IF ( iom_use( "ymtrpsnw" ) ) CALL iom_put( "ymtrpsnw" , diag_ymtrp_snw(:,:) ) ! Y-component of snow mass transport346 347 IF ( iom_use( "xatrp" ) ) CALL iom_put( "xatrp" , diag_xatrp(:,:) ) ! X-component of ice area transport348 IF ( iom_use( "yatrp" ) ) CALL iom_put( "yatrp" , diag_yatrp(:,:) ) ! Y-component of ice area transport349 350 305 IF ( iom_use( "utau_ice" ) ) CALL iom_put( "utau_ice" , utau_ice(:,:) * zswi(:,:) + zmiss(:,:) ) ! Wind stress term in force balance (x) 351 306 IF ( iom_use( "vtau_ice" ) ) CALL iom_put( "vtau_ice" , vtau_ice(:,:) * zswi(:,:) + zmiss(:,:) ) ! Wind stress term in force balance (y) 352 307 353 IF ( iom_use( "utau_oi" ) ) CALL iom_put( "utau_oi" , diag_utau_oi(:,:) * zswi(:,:) + zmiss(:,:) ) ! Ocean stress term in force balance (x)354 IF ( iom_use( "vtau_oi" ) ) CALL iom_put( "vtau_oi" , diag_vtau_oi(:,:) * zswi(:,:) + zmiss(:,:) ) ! Ocean stress term in force balance (y)355 308 356 309 IF ( iom_use( "icestr" ) ) CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) + zmiss(:,:) ) ! Ice strength 357 358 IF ( iom_use( "dssh_dx" ) ) CALL iom_put( "dssh_dx" , diag_dssh_dx(:,:) * zswi(:,:) + zmiss(:,:) ) ! Sea-surface tilt term in force balance (x)359 IF ( iom_use( "dssh_dy" ) ) CALL iom_put( "dssh_dy" , diag_dssh_dy(:,:) * zswi(:,:) + zmiss(:,:) ) ! Sea-surface tilt term in force balance (y)360 361 IF ( iom_use( "corstrx" ) ) CALL iom_put( "corstrx" , diag_corstrx(:,:) * zswi(:,:) + zmiss(:,:) ) ! Coriolis force term in force balance (x)362 IF ( iom_use( "corstry" ) ) CALL iom_put( "corstry" , diag_corstry(:,:) * zswi(:,:) + zmiss(:,:) ) ! Coriolis force term in force balance (y)363 364 IF ( iom_use( "intstrx" ) ) CALL iom_put( "intstrx" , diag_intstrx(:,:) * zswi(:,:) + zmiss(:,:) ) ! Internal force term in force balance (x)365 IF ( iom_use( "intstry" ) ) CALL iom_put( "intstry" , diag_intstry(:,:) * zswi(:,:) + zmiss(:,:) ) ! Internal force term in force balance (y)366 367 IF ( iom_use( "normstr" ) ) CALL iom_put( "normstr" , diag_sig1(:,:) * zswi(:,:) + zmiss(:,:) ) ! Normal stress368 IF ( iom_use( "sheastr" ) ) CALL iom_put( "sheastr" , diag_sig2(:,:) * zswi(:,:) + zmiss(:,:) ) ! Shear stress369 310 370 311 !--------------------------------
Note: See TracChangeset
for help on using the changeset viewer.