diff --git a/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg b/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg index 8e1b5a14cbd2cc8341fbf8c982e5f71f30495457..31af4cd1171848e621f1a1ee60371649f3298ac2 100644 --- a/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg +++ b/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg @@ -325,7 +325,7 @@ sn_sal = 'dyna_grid_T' , 120. , 'vosaline' , .true. , .true. , 'yearly' , '' , '' , '' sn_mld = 'dyna_grid_T' , 120. , 'somixhgt' , .true. , .true. , 'yearly' , '' , '' , '' sn_emp = 'dyna_grid_T' , 120. , 'sowaflcd' , .true. , .true. , 'yearly' , '' , '' , '' - sn_fmf = 'dyna_grid_T' , 120. , 'iowaflup' , .true. , .true. , 'yearly' , '' , '' , '' + sn_fwf = 'dyna_grid_T' , 120. , 'iowaflup' , .true. , .true. , 'yearly' , '' , '' , '' sn_ice = 'dyna_grid_T' , 120. , 'soicecov' , .true. , .true. , 'yearly' , '' , '' , '' sn_qsr = 'dyna_grid_T' , 120. , 'soshfldo' , .true. , .true. , 'yearly' , '' , '' , '' sn_wnd = 'dyna_grid_T' , 120. , 'sowindsp' , .true. , .true. , 'yearly' , '' , '' , '' diff --git a/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg b/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg index e365edaa458129637a1d1f2a466f7cebc526a0b4..380fbc570eb6f76dce449b9fce4cb66f630a51ab 100644 --- a/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg +++ b/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg @@ -324,7 +324,7 @@ sn_mld = 'dyna_grid_T' , 120. , 'mldr10_1' , .true. , .true. , 'yearly' , '' , '' , '' sn_emp = 'dyna_grid_T' , 120. , 'wfo' , .true. , .true. , 'yearly' , '' , '' , '' sn_empb = 'dyna_grid_T' , 120. , 'wfob' , .true. , .true. , 'yearly' , '' , '' , '' - sn_fmf = 'dyna_grid_T' , 120. , 'fmmflx' , .true. , .true. , 'yearly' , '' , '' , '' + sn_fwf = 'dyna_grid_T' , 120. , 'iowaflup' , .true. , .true. , 'yearly' , '' , '' , '' sn_rnf = 'dyna_grid_T' , 120. , 'runoffs' , .true. , .true. , 'yearly' , '' , '' , '' sn_ice = 'dyna_grid_T' , 120. , 'siconc' , .true. , .true. , 'yearly' , '' , '' , '' sn_qsr = 'dyna_grid_T' , 120. , 'rsntds' , .true. , .true. , 'yearly' , '' , '' , '' diff --git a/cfgs/SHARED/field_def_nemo-oce.xml b/cfgs/SHARED/field_def_nemo-oce.xml index 8829cd382a084db3361a17e39d0ee4337cfc59ac..6f71577972d49e930714b1cf5ee0e79d445880c5 100644 --- a/cfgs/SHARED/field_def_nemo-oce.xml +++ b/cfgs/SHARED/field_def_nemo-oce.xml @@ -415,7 +415,7 @@ that are available in the tidal-forcing implementation (see - + diff --git a/cfgs/SHARED/namelist_ref b/cfgs/SHARED/namelist_ref index 8cca2df5f1883b2d4cf966d269383d67068a1249..c235739131b153abb3d43dc918033436878c0de6 100644 --- a/cfgs/SHARED/namelist_ref +++ b/cfgs/SHARED/namelist_ref @@ -1117,7 +1117,7 @@ sn_sal = 'dyna_grid_T' , 120. , 'vosaline' , .true. , .true. , 'yearly' , '' , '' , '' sn_mld = 'dyna_grid_T' , 120. , 'somixhgt' , .true. , .true. , 'yearly' , '' , '' , '' sn_emp = 'dyna_grid_T' , 120. , 'sowaflup' , .true. , .true. , 'yearly' , '' , '' , '' - sn_fmf = 'dyna_grid_T' , 120. , 'iowaflup' , .true. , .true. , 'yearly' , '' , '' , '' + sn_fwf = 'dyna_grid_T' , 120. , 'iowaflup' , .true. , .true. , 'yearly' , '' , '' , '' sn_ice = 'dyna_grid_T' , 120. , 'soicecov' , .true. , .true. , 'yearly' , '' , '' , '' sn_qsr = 'dyna_grid_T' , 120. , 'soshfldo' , .true. , .true. , 'yearly' , '' , '' , '' sn_wnd = 'dyna_grid_T' , 120. , 'sowindsp' , .true. , .true. , 'yearly' , '' , '' , '' diff --git a/src/ICE/ice.F90 b/src/ICE/ice.F90 index 8e5e48747f87453c754d850c588b21b8591e7f0a..1eb9d9943214f7a294fb9381f440d996ee0b059b 100644 --- a/src/ICE/ice.F90 +++ b/src/ICE/ice.F90 @@ -257,7 +257,6 @@ MODULE ice REAL(wp), PUBLIC :: r1_Dt_ice !: = 1. / rDt_ice REAL(wp), PUBLIC :: r1_nlay_i !: 1 / nlay_i REAL(wp), PUBLIC :: r1_nlay_s !: 1 / nlay_s - REAL(wp), PUBLIC :: rswitch !: switch for the presence of ice (1) or not (0) REAL(wp), PUBLIC :: rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft !: conservation diagnostics REAL(wp), PUBLIC, PARAMETER :: epsi06 = 1.e-06_wp !: small number REAL(wp), PUBLIC, PARAMETER :: epsi10 = 1.e-10_wp !: small number diff --git a/src/ICE/icedyn_adv_pra.F90 b/src/ICE/icedyn_adv_pra.F90 index 5bbde365e7f22c3c53ea3abfa2f28f9253d743dd..fff6785e4a12977b4fe16f1737904f23fcb4be69 100644 --- a/src/ICE/icedyn_adv_pra.F90 +++ b/src/ICE/icedyn_adv_pra.F90 @@ -449,7 +449,6 @@ CONTAINS REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - REAL(wp) :: zpsm, zps0 REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy - REAL(wp) :: zswitch REAL(wp), DIMENSION(jpi,jpj) :: zf0 , zfx , zfy , zbet ! 2D workspace REAL(wp), DIMENSION(jpi,jpj) :: zfm , zfxx , zfyy , zfxy ! - - REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - @@ -472,32 +471,40 @@ CONTAINS zpsxy = psxy(ji,jj,jcat) ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) - zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) + zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1._wp - pcrh ) * zpsm , epsi20 ) ! zslpmax = MAX( 0._wp, zps0 ) - zs1max = 1.5 * zslpmax + zs1max = 1.5_wp * zslpmax zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) ) - zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) - zswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask - + zs2new = MIN( 2._wp * zslpmax - 0.3334_wp * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) + ! zps0 = zslpmax - zpsx = zs1new * zswitch - zpsxx = zs2new * zswitch - zpsy = zpsy * zswitch - zpsyy = zpsyy * zswitch - zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * zswitch - + ! + IF( zslpmax > 0._wp ) THEN + zpsx = zs1new * tmask(ji,jj,1) + zpsxx = zs2new * tmask(ji,jj,1) + zpsy = zpsy * tmask(ji,jj,1) + zpsyy = zpsyy * tmask(ji,jj,1) + zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * tmask(ji,jj,1) + ELSE + zpsx = 0._wp + zpsxx = 0._wp + zpsy = 0._wp + zpsyy = 0._wp + zpsxy = 0._wp + ENDIF + ! ! Calculate fluxes and moments between boxes i<-->i+1 ! ! Flux from i to i+1 WHEN u GT 0 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm zalfq = zalf * zalf - zalf1 = 1.0 - zalf + zalf1 = 1._wp - zalf zalf1q = zalf1 * zalf1 ! zfm (ji,jj) = zalf * zpsm zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) - zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx ) + zfx (ji,jj) = zalfq * ( zpsx + 3._wp * zalf1 * zpsxx ) zfxx(ji,jj) = zalf * zpsxx * zalfq zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy ) zfxy(ji,jj) = zalfq * zpsxy @@ -506,7 +513,7 @@ CONTAINS ! ! Readjust moments remaining in the box. zpsm = zpsm - zfm(ji,jj) zps0 = zps0 - zf0(ji,jj) - zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) + zpsx = zalf1q * ( zpsx - 3._wp * zalf * zpsxx ) zpsxx = zalf1 * zalf1q * zpsxx zpsy = zpsy - zfy (ji,jj) zpsyy = zpsyy - zfyy(ji,jj) @@ -527,7 +534,7 @@ CONTAINS zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj) zalg (ji,jj) = zalf zalfq = zalf * zalf - zalf1 = 1.0 - zalf + zalf1 = 1._wp - zalf zalg1 (ji,jj) = zalf1 zalf1q = zalf1 * zalf1 zalg1q(ji,jj) = zalf1q @@ -535,7 +542,7 @@ CONTAINS zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj) zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj) & & - zalf1 * ( psx(ji+1,jj,jcat) - (zalf1 - zalf ) * psxx(ji+1,jj,jcat) ) ) - zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jcat) - 3.0 * zalf1 * psxx(ji+1,jj,jcat) ) + zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jcat) - 3._wp * zalf1 * psxx(ji+1,jj,jcat) ) zfxx (ji,jj) = zfxx(ji,jj) + zalf * ( psxx(ji+1,jj,jcat) * zalfq ) zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jcat) - zalf1 * psxy(ji+1,jj,jcat) ) zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jcat) @@ -552,12 +559,12 @@ CONTAINS zpsyy = psyy(ji,jj,jcat) zpsxy = psxy(ji,jj,jcat) ! ! Readjust moments remaining in the box. - zbt = zbet(ji-1,jj) - zbt1 = 1.0 - zbet(ji-1,jj) + zbt = zbet(ji-1,jj) + zbt1 = 1._wp - zbet(ji-1,jj) ! zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) - zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) + zpsx = zalg1q(ji-1,jj) * ( zpsx + 3._wp * zalg(ji-1,jj) * zpsxx ) zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) ) zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) @@ -565,38 +572,38 @@ CONTAINS ! Put the temporary moments into appropriate neighboring boxes. ! ! Flux from i to i+1 IF u GT 0. - zbt = zbet(ji-1,jj) - zbt1 = 1.0 - zbet(ji-1,jj) + zbt = zbet(ji-1,jj) + zbt1 = 1._wp - zbet(ji-1,jj) zpsm = zbt1 * zpsm + zbt * ( zpsm + zfm(ji-1,jj) ) zalf = zbt * zfm(ji-1,jj) / zpsm - zalf1 = 1.0 - zalf + zalf1 = 1._wp - zalf ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) ! zps0 = zbt1 * zps0 + zbt * ( zps0 + zf0(ji-1,jj) ) - zpsx = zbt1 * zpsx + zbt * ( ( zalf * zfx(ji-1,jj) + zalf1 * zpsx ) + 3.0 * ztemp ) + zpsx = zbt1 * zpsx + zbt * ( ( zalf * zfx(ji-1,jj) + zalf1 * zpsx ) + 3._wp * ztemp ) zpsxx = zbt1 * zpsxx + zbt * ( ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx ) & - & + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) + & + 5._wp * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) zpsxy = zbt1 * zpsxy + zbt * ( ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy ) & - & + 3.0 * ( -zalf1 * zfy(ji-1,jj) + zalf * zpsy ) ) + & + 3._wp * ( -zalf1 * zfy(ji-1,jj) + zalf * zpsy ) ) zpsy = zbt1 * zpsy + zbt * ( zpsy + zfy (ji-1,jj) ) zpsyy = zbt1 * zpsyy + zbt * ( zpsyy + zfyy(ji-1,jj) ) ! ! Flux from i+1 to i IF u LT 0. - zbt = zbet(ji,jj) - zbt1 = 1.0 - zbet(ji,jj) - zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji ,jj) ) - zalf = zbt1 * zfm(ji ,jj) / zpsm - zalf1 = 1.0 - zalf - ztemp = -zalf * zps0 + zalf1 * zf0(ji ,jj) + zbt = zbet(ji,jj) + zbt1 = 1._wp - zbet(ji,jj) + zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) + zalf = zbt1 * zfm(ji,jj) / zpsm + zalf1 = 1._wp - zalf + ztemp = -zalf * zps0 + zalf1 * zf0(ji,jj) ! - zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji ,jj) ) - zpsx = zbt * zpsx + zbt1 * ( ( zalf * zfx(ji ,jj) + zalf1 * zpsx ) + 3.0 * ztemp ) - zpsxx = zbt * zpsxx + zbt1 * ( ( zalf * zalf * zfxx(ji ,jj) + zalf1 * zalf1 * zpsxx ) & - & + 5.0 * ( zalf * zalf1 * ( -zpsx + zfx(ji ,jj) ) + ( zalf1 - zalf ) * ztemp ) ) - zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji ,jj) + zalf1 * zpsxy ) & - & + 3.0 * ( zalf1 * zfy(ji ,jj) - zalf * zpsy ) ) - zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji ,jj) ) - zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji ,jj) ) + zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) ) + zpsx = zbt * zpsx + zbt1 * ( ( zalf * zfx(ji,jj) + zalf1 * zpsx ) + 3._wp * ztemp ) + zpsxx = zbt * zpsxx + zbt1 * ( ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx ) & + & + 5._wp * ( zalf * zalf1 * ( -zpsx + zfx(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) + zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji,jj) + zalf1 * zpsxy ) & + & + 3._wp * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) + zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji,jj) ) + zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) ! psm (ji,jj) = zpsm ! optimization ps0 (ji,jj) = zps0 @@ -636,7 +643,6 @@ CONTAINS REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - REAL(wp) :: zpsm, zps0 REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy - REAL(wp) :: zswitch REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - - REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - @@ -659,32 +665,40 @@ CONTAINS zpsxy = psxy(ji,jj,jcat) ! ! Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) - zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) + zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1._wp - pcrh ) * zpsm , epsi20 ) ! zslpmax = MAX( 0._wp, zps0 ) - zs1max = 1.5 * zslpmax + zs1max = 1.5_wp * zslpmax zs1new = MIN( zs1max, MAX( -zs1max, zpsy ) ) - zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new ) - zslpmax, zpsyy ) ) - zswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask + zs2new = MIN( ( 2._wp * zslpmax - 0.3334_wp * ABS( zs1new ) ), MAX( ABS( zs1new ) - zslpmax, zpsyy ) ) ! zps0 = zslpmax - zpsx = zpsx * zswitch - zpsxx = zpsxx * zswitch - zpsy = zs1new * zswitch - zpsyy = zs2new * zswitch - zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * zswitch - + ! + IF( zslpmax > 0._wp ) THEN + zpsx = zpsx * tmask(ji,jj,1) + zpsxx = zpsxx * tmask(ji,jj,1) + zpsy = zs1new * tmask(ji,jj,1) + zpsyy = zs2new * tmask(ji,jj,1) + zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * tmask(ji,jj,1) + ELSE + zpsx = 0._wp + zpsxx = 0._wp + zpsy = 0._wp + zpsyy = 0._wp + zpsxy = 0._wp + ENDIF + ! ! Calculate fluxes and moments between boxes j<-->j+1 ! ! Flux from j to j+1 WHEN v GT 0 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm zalfq = zalf * zalf - zalf1 = 1.0 - zalf + zalf1 = 1._wp - zalf zalf1q = zalf1 * zalf1 ! zfm (ji,jj) = zalf * zpsm zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1 - zalf) * zpsyy ) ) - zfy (ji,jj) = zalfq * ( zpsy + 3.0 * zalf1 * zpsyy ) + zfy (ji,jj) = zalfq * ( zpsy + 3._wp * zalf1 * zpsyy ) zfyy(ji,jj) = zalf * zpsyy * zalfq zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy ) zfxy(ji,jj) = zalfq * zpsxy @@ -693,7 +707,7 @@ CONTAINS ! ! Readjust moments remaining in the box. zpsm = zpsm - zfm(ji,jj) zps0 = zps0 - zf0(ji,jj) - zpsy = zalf1q * ( zpsy -3.0 * zalf * zpsyy ) + zpsy = zalf1q * ( zpsy -3._wp * zalf * zpsyy ) zpsyy = zalf1 * zalf1q * zpsyy zpsx = zpsx - zfx(ji,jj) zpsxx = zpsxx - zfxx(ji,jj) @@ -713,7 +727,7 @@ CONTAINS zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1) zalg (ji,jj) = zalf zalfq = zalf * zalf - zalf1 = 1.0 - zalf + zalf1 = 1._wp - zalf zalg1 (ji,jj) = zalf1 zalf1q = zalf1 * zalf1 zalg1q(ji,jj) = zalf1q @@ -721,7 +735,7 @@ CONTAINS zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1) zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1) & & - zalf1 * (psy(ji,jj+1,jcat) - (zalf1 - zalf ) * psyy(ji,jj+1,jcat) ) ) - zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jcat) - 3.0 * zalf1 * psyy(ji,jj+1,jcat) ) + zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jcat) - 3._wp * zalf1 * psyy(ji,jj+1,jcat) ) zfyy (ji,jj) = zfyy(ji,jj) + zalf * ( psyy(ji,jj+1,jcat) * zalfq ) zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jcat) - zalf1 * psxy(ji,jj+1,jcat) ) zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jcat) @@ -739,11 +753,11 @@ CONTAINS zpsxy = psxy(ji,jj,jcat) ! ! Readjust moments remaining in the box. zbt = zbet(ji,jj-1) - zbt1 = ( 1.0 - zbet(ji,jj-1) ) + zbt1 = ( 1._wp - zbet(ji,jj-1) ) ! zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) - zpsy = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) + zpsy = zalg1q(ji,jj-1) * ( zpsy + 3._wp * zalg(ji,jj-1) * zpsyy ) zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) ) zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) @@ -752,37 +766,37 @@ CONTAINS ! Put the temporary moments into appropriate neighboring boxes. ! ! Flux from j to j+1 IF v GT 0. zbt = zbet(ji,jj-1) - zbt1 = 1.0 - zbet(ji,jj-1) + zbt1 = 1._wp - zbet(ji,jj-1) zpsm = zbt1 * zpsm + zbt * ( zpsm + zfm(ji,jj-1) ) zalf = zbt * zfm(ji,jj-1) / zpsm - zalf1 = 1.0 - zalf + zalf1 = 1._wp - zalf ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) ! zps0 = zbt1 * zps0 + zbt * ( zps0 + zf0(ji,jj-1) ) - zpsy = zbt1 * zpsy + zbt * ( ( zalf * zfy(ji,jj-1) + zalf1 * zpsy ) + 3.0 * ztemp ) + zpsy = zbt1 * zpsy + zbt * ( ( zalf * zfy(ji,jj-1) + zalf1 * zpsy ) + 3._wp * ztemp ) zpsyy = zbt1 * zpsyy + zbt * ( ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy ) & - & + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) + & + 5._wp * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) zpsxy = zbt1 * zpsxy + zbt * ( ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy ) & - & + 3.0 * ( -zalf1 * zfx(ji,jj-1) + zalf * zpsx ) ) + & + 3._wp * ( -zalf1 * zfx(ji,jj-1) + zalf * zpsx ) ) zpsx = zbt1 * zpsx + zbt * ( zpsx + zfx (ji,jj-1) ) zpsxx = zbt1 * zpsxx + zbt * ( zpsxx + zfxx(ji,jj-1) ) ! ! Flux from j+1 to j IF v LT 0. zbt = zbet(ji,jj) - zbt1 = 1.0 - zbet(ji,jj) + zbt1 = 1._wp - zbet(ji,jj) zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) zalf = zbt1 * zfm(ji,jj) / zpsm - zalf1 = 1.0 - zalf + zalf1 = 1._wp - zalf ztemp = -zalf * zps0 + zalf1 * zf0(ji,jj) ! - zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj ) ) - zpsy = zbt * zpsy + zbt1 * ( ( zalf * zfy(ji,jj ) + zalf1 * zpsy ) + 3.0 * ztemp ) - zpsyy = zbt * zpsyy + zbt1 * ( ( zalf * zalf * zfyy(ji,jj ) + zalf1 * zalf1 * zpsyy ) & - & + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj ) ) + ( zalf1 - zalf ) * ztemp ) ) - zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji,jj ) + zalf1 * zpsxy ) & - & + 3.0 * ( zalf1 * zfx(ji,jj ) - zalf * zpsx ) ) - zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj ) ) - zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj ) ) + zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) ) + zpsy = zbt * zpsy + zbt1 * ( ( zalf * zfy(ji,jj) + zalf1 * zpsy ) + 3._wp * ztemp ) + zpsyy = zbt * zpsyy + zbt1 * ( ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy ) & + & + 5._wp * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) + zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji,jj) + zalf1 * zpsxy ) & + & + 3._wp * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) + zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj) ) + zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) ! psm (ji,jj) = zpsm ! optimization ps0 (ji,jj) = zps0 diff --git a/src/ICE/icedyn_rdgrft.F90 b/src/ICE/icedyn_rdgrft.F90 index 33fe074d63528bc212d2fef7443f0a55fa94bc0a..ad7d1de81c22ea810c7dc50eb8257aa2b7da1db2 100644 --- a/src/ICE/icedyn_rdgrft.F90 +++ b/src/ICE/icedyn_rdgrft.F90 @@ -839,7 +839,6 @@ CONTAINS INTEGER :: ji, jj, jl ! dummy loop indices REAL(wp) :: z1_3 ! local scalars REAL(wp), DIMENSION(A2D(0)) :: zworka ! temporary array used here - REAL(wp), DIMENSION(jpi,jpj) :: zmsk ! temporary array used here !! LOGICAL :: ln_str_R75 REAL(wp) :: zhi, zcp @@ -847,11 +846,8 @@ CONTAINS REAL(wp), PARAMETER :: zmax_strength = 200.e3_wp ! Richter-Menge and Elder (1998) estimate maximum in Beaufort Sea in wintertime of the order 150 kN/m REAL(wp), DIMENSION(jpij) :: zstrength, zaksum ! strength in 1D !!---------------------------------------------------------------------- - ! prepare the mask + ! at_i needed for strength at_i(:,:) = SUM( a_i, dim=3 ) - DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) - zmsk(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 1 if ice , 0 if no ice - END_2D ! SELECT CASE( nice_str ) !--- Set which ice strength is chosen @@ -940,20 +936,28 @@ CONTAINS CASE ( np_strh79 ) !== Hibler(1979)'s method ==! ! DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) - strength(ji,jj) = rn_pstar * SUM( v_i(ji,jj,:) ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) * zmsk(ji,jj) + IF( at_i(ji,jj) > epsi10 ) THEN + strength(ji,jj) = rn_pstar * SUM( v_i(ji,jj,:) ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) + ELSE + strength(ji,jj) = 0._wp + ENDIF END_2D ! CASE ( np_strcst ) !== Constant strength ==! ! DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) - strength(ji,jj) = rn_str * zmsk(ji,jj) + IF( at_i(ji,jj) > epsi10 ) THEN + strength(ji,jj) = rn_str + ELSE + strength(ji,jj) = 0._wp + ENDIF END_2D ! END SELECT ! IF( ln_str_smooth ) THEN !--- Spatial smoothing DO_2D( 0, 0, 0, 0 ) - IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN + IF( at_i(ji,jj) > epsi10 ) THEN zworka(ji,jj) = ( 4._wp * strength(ji,jj) & & + ( ( strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) ) & & + ( strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) ) ) & @@ -964,7 +968,7 @@ CONTAINS END_2D DO_2D( 0, 0, 0, 0 ) - strength(ji,jj) = zworka(ji,jj) * zmsk(ji,jj) + strength(ji,jj) = zworka(ji,jj) END_2D CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp ) ! diff --git a/src/ICE/icedyn_rhg_evp.F90 b/src/ICE/icedyn_rhg_evp.F90 index 59f9ef866e95f427f46f285553ca0fec8eef62c5..46898cf6b8f307dca435c2832a99035f3c50bcdf 100644 --- a/src/ICE/icedyn_rhg_evp.F90 +++ b/src/ICE/icedyn_rhg_evp.F90 @@ -135,8 +135,6 @@ CONTAINS ! REAL(wp) :: zfac_x, zfac_y ! - REAL(wp) :: zswitch - ! REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta, P/delta at T points REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 ! @@ -190,15 +188,18 @@ CONTAINS IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' ! DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) - zmsk (ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 1 if ice , 0 if no ice + IF( at_i(ji,jj) < epsi10 ) THEN ; zmsk(ji,jj) = 0._wp + ELSE ; zmsk(ji,jj) = 1._wp ; ENDIF END_2D ! for diagnostics and convergence tests DO_2D( 0, 0, 0, 0 ) - zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice + IF( at_i(ji,jj) < epsi06 ) THEN ; zmsk00(ji,jj) = 0._wp + ELSE ; zmsk00(ji,jj) = 1._wp ; ENDIF END_2D IF( nn_rhg_chkcvg > 0 ) THEN DO_2D( 0, 0, 0, 0 ) - zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less + IF( at_i(ji,jj) < 0.15_wp ) THEN ; zmsk15(ji,jj) = 0._wp + ELSE ; zmsk15(ji,jj) = 1._wp ; ENDIF END_2D ENDIF ! @@ -313,9 +314,11 @@ CONTAINS zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) ! masks - zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice - zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice - + IF( zmassU > 0._wp ) THEN ; zmsk00x(ji,jj) = 1._wp + ELSE ; zmsk00x(ji,jj) = 0._wp ; ENDIF + IF( zmassV > 0._wp ) THEN ; zmsk00y(ji,jj) = 1._wp + ELSE ; zmsk00y(ji,jj) = 0._wp ; ENDIF + ! switches IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF @@ -472,17 +475,17 @@ CONTAINS ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility) DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! !--- U points - zfU(ji,jj) = 0.5_wp * ( ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & - & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & - & ) * r1_e2u(ji,jj) ) & + zfU(ji,jj) = 0.5_wp * ( ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & + & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & + & ) * r1_e2u(ji,jj) ) & & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & & ) * 2._wp * r1_e1u(ji,jj) & & ) * r1_e1e2u(ji,jj) ! ! !--- V points - zfV(ji,jj) = 0.5_wp * ( ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & - & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & - & ) * r1_e1v(ji,jj) ) & + zfV(ji,jj) = 0.5_wp * ( ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & + & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & + & ) * r1_e1v(ji,jj) ) & & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & & ) * 2._wp * r1_e2v(ji,jj) & & ) * r1_e1e2v(ji,jj) @@ -496,6 +499,15 @@ CONTAINS ! --- Computation of ice velocity --- ! ! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017 ! Bouillon et al. 2009 (eq 34-35) => stable + ! + ! aEVP formulation given by eq.6 from Kimmritz et al. 2016, but taken the last term implicitly (as in eq. 8) + ! u(p+1) - u(p) = 1/beta * ( dt/m * RHS(p+1) + u(n) - u(p+1) ) + ! with RHS = tau_ai + tau_oi + tau_bi + Coriolis + Spg + ! and tau_oi = ztau0 * ( uoce - u(p+1) ) + ! tau_bi = ztauB * u(p+1) + ! Hence: + ! u(p+1) = 1/(m/dt*(beta+1)+ztauO-ztauB) * (m/dt*(beta*u(p)+u(n))+RHS+ztauO*u(p)) + ! IF( MOD(jter,2) == 0 ) THEN ! even iterations ! DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) @@ -519,29 +531,32 @@ CONTAINS ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) ! - ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) - ! 1 = sliding friction : TauB < RHS - zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) - ! - IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) + IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) + ! zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) - v_ice(ji,jj) = ( ( zswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity - & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) - & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast - & + ( 1._wp - zswitch ) * ( v_ice_b(ji,jj) & - & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 - & ) / ( zbetav + 1._wp ) & - & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin - & ) * zmsk00y(ji,jj) - ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) - v_ice(ji,jj) = ( ( zswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity - & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) - & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast - & + ( 1._wp - zswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 - & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin - & ) * zmsk00y(ji,jj) + ! + IF( ( zRHS + ztauy_base(ji,jj) ) < 0._wp .AND. zRHS >= 0._wp ) THEN ! static friction => slow decrease to v=0 + v_ice(ji,jj) = ( v_ice_b(ji,jj) + v_ice(ji,jj) * MAX( 0._wp, zbetav - zdtevp*rn_lf_relax ) ) / (zbetav+1._wp) + ELSE + v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) + zRHS + zTauO * v_ice(ji,jj) ) & + & / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) + ENDIF + ! + ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) + ! + IF( ( zRHS + ztauy_base(ji,jj) ) < 0._wp .AND. zRHS >= 0._wp ) THEN ! static friction => slow decrease to v=0 + v_ice(ji,jj) = v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) + ELSE + v_ice(ji,jj) = ( zmV_t(ji,jj) * v_ice(ji,jj) + zRHS + zTauO * v_ice(ji,jj) ) & + & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) + ENDIF + ! ENDIF + ! v_ice = v_oce/100 if mass < zmmin & conc < zamin + v_ice(ji,jj) = ( v_ice(ji,jj)*zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * (1._wp - zmsk01y(ji,jj)) ) * zmsk00y(ji,jj) + ! END_2D + ! IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) ! DO_2D( 0, 0, 0, 0 ) @@ -565,29 +580,32 @@ CONTAINS ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) ! - ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) - ! 1 = sliding friction : TauB < RHS - zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) - ! IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) + ! zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) - u_ice(ji,jj) = ( ( zswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity - & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) - & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast - & + ( 1._wp - zswitch ) * ( u_ice_b(ji,jj) & - & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 - & ) / ( zbetau + 1._wp ) & - & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin - & ) * zmsk00x(ji,jj) + ! + IF( ( zRHS + ztaux_base(ji,jj) ) < 0._wp .AND. zRHS >= 0._wp ) THEN ! static friction => slow decrease to v=0 + u_ice(ji,jj) = ( u_ice_b(ji,jj) + u_ice(ji,jj) * MAX( 0._wp, zbetau - zdtevp*rn_lf_relax ) ) / (zbetau+1._wp) + ELSE + u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) + zRHS + zTauO * u_ice(ji,jj) ) & + & / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) + ENDIF + ! ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) - u_ice(ji,jj) = ( ( zswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity - & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) - & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast - & + ( 1._wp - zswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 - & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin - & ) * zmsk00x(ji,jj) + ! + IF( ( zRHS + ztaux_base(ji,jj) ) < 0._wp .AND. zRHS >= 0._wp ) THEN ! static friction => slow decrease to v=0 + u_ice(ji,jj) = u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) + ELSE + u_ice(ji,jj) = ( zmU_t(ji,jj) * u_ice(ji,jj) + zRHS + zTauO * u_ice(ji,jj) ) & + & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) + ENDIF + ! ENDIF + ! u_ice = u_oce/100 if mass < zmmin & conc < zamin + u_ice(ji,jj) = ( u_ice(ji,jj)*zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * (1._wp - zmsk01x(ji,jj)) ) * zmsk00x(ji,jj) + ! END_2D + ! IF( nn_hls == 1 ) THEN ; CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) ELSE ; CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) ENDIF @@ -615,29 +633,32 @@ CONTAINS ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) ! - ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) - ! 1 = sliding friction : TauB < RHS - zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) - ! IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) + ! zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) - u_ice(ji,jj) = ( ( zswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity - & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) - & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast - & + ( 1._wp - zswitch ) * ( u_ice_b(ji,jj) & - & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 - & ) / ( zbetau + 1._wp ) & - & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin - & ) * zmsk00x(ji,jj) + ! + IF( ( zRHS + ztaux_base(ji,jj) ) < 0._wp .AND. zRHS >= 0._wp ) THEN ! static friction => slow decrease to v=0 + u_ice(ji,jj) = ( u_ice_b(ji,jj) + u_ice(ji,jj) * MAX( 0._wp, zbetau - zdtevp*rn_lf_relax ) ) / (zbetau+1._wp) + ELSE + u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) + zRHS + zTauO * u_ice(ji,jj) ) & + & / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) + ENDIF + ! ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) - u_ice(ji,jj) = ( ( zswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity - & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) - & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast - & + ( 1._wp - zswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 - & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin - & ) * zmsk00x(ji,jj) + ! + IF( ( zRHS + ztaux_base(ji,jj) ) < 0._wp .AND. zRHS >= 0._wp ) THEN ! static friction => slow decrease to v=0 + u_ice(ji,jj) = u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) + ELSE + u_ice(ji,jj) = ( zmU_t(ji,jj) * u_ice(ji,jj) + zRHS + zTauO * u_ice(ji,jj) ) & + & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) + ENDIF + ! ENDIF + ! u_ice = u_oce/100 if mass < zmmin & conc < zamin + u_ice(ji,jj) = ( u_ice(ji,jj)*zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * (1._wp - zmsk01x(ji,jj)) ) * zmsk00x(ji,jj) + ! END_2D + ! IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) ! DO_2D( 0, 0, 0, 0 ) @@ -661,29 +682,32 @@ CONTAINS ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) ! - ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) - ! 1 = sliding friction : TauB < RHS - zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) - ! IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) + ! zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) - v_ice(ji,jj) = ( ( zswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity - & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) - & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast - & + ( 1._wp - zswitch ) * ( v_ice_b(ji,jj) & - & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 - & ) / ( zbetav + 1._wp ) & - & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin - & ) * zmsk00y(ji,jj) + ! + IF( ( zRHS + ztauy_base(ji,jj) ) < 0._wp .AND. zRHS >= 0._wp ) THEN ! static friction => slow decrease to v=0 + v_ice(ji,jj) = ( v_ice_b(ji,jj) + v_ice(ji,jj) * MAX( 0._wp, zbetav - zdtevp*rn_lf_relax ) ) / (zbetav+1._wp) + ELSE + v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) + zRHS + zTauO * v_ice(ji,jj) ) & + & / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) + ENDIF + ! ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) - v_ice(ji,jj) = ( ( zswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity - & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) - & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast - & + ( 1._wp - zswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 - & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin - & ) * zmsk00y(ji,jj) + ! + IF( ( zRHS + ztauy_base(ji,jj) ) < 0._wp .AND. zRHS >= 0._wp ) THEN ! static friction => slow decrease to v=0 + v_ice(ji,jj) = v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) + ELSE + v_ice(ji,jj) = ( zmV_t(ji,jj) * v_ice(ji,jj) + zRHS + zTauO * v_ice(ji,jj) ) & + & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) + ENDIF + ! ENDIF + ! v_ice = v_oce/100 if mass < zmmin & conc < zamin + v_ice(ji,jj) = ( v_ice(ji,jj)*zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * (1._wp - zmsk01y(ji,jj)) ) * zmsk00y(ji,jj) + ! END_2D + ! IF( nn_hls == 1 ) THEN ; CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) ELSE ; CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) ENDIF @@ -775,12 +799,12 @@ CONTAINS zdelta(ji,jj) = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! delta ! delta* at T points (pdelta_i) - zswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0 - pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * zswitch + IF( zdelta(ji,jj) > 0._wp ) THEN ; pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl + ELSE ; pdelta_i(ji,jj) = 0._wp + ENDIF ! it seems that deformation used for advection and mech redistribution is delta* ! MV in principle adding creep limit is a regularization for viscosity not for delta ! delta_star should not (in my view) be used in a replacement for delta - END_2D CALL lbc_lnk( 'icedyn_rhg_evp', zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp ) diff --git a/src/ICE/icedyn_rhg_vp.F90 b/src/ICE/icedyn_rhg_vp.F90 index 1553424ea429d9093bdefb34aeca86ab0887dddc..53073a428a35e5f8dd37e56e7b26db449c3353ff 100644 --- a/src/ICE/icedyn_rhg_vp.F90 +++ b/src/ICE/icedyn_rhg_vp.F90 @@ -142,6 +142,7 @@ CONTAINS INTEGER :: nn_zebra_vp ! number of zebra steps ! + REAL(wp) :: zswitch REAL(wp) :: zrhoco ! rho0 * rn_cio REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity REAL(wp) :: zglob_area ! global ice area for diagnostics @@ -1097,8 +1098,8 @@ CONTAINS zdelta(ji,jj) = zfac ! delta* at T points - rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 - pdelta_i(ji,jj) = zfac + rn_creepl ! * rswitch + zswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 + pdelta_i(ji,jj) = zfac + rn_creepl ! * zswitch END_2D diff --git a/src/ICE/icesbc.F90 b/src/ICE/icesbc.F90 index 2b94b08208cf667241ecfc6a8e490a5af3f25845..2e8baf6bfd57958bd25494734351129a59c86954 100644 --- a/src/ICE/icesbc.F90 +++ b/src/ICE/icesbc.F90 @@ -324,7 +324,9 @@ CONTAINS ! Partial computation of forcing for the thermodynamic sea ice model !--------------------------------------------------------------------! DO_2D( 0, 0, 0, 0 ) ! needed for qlead - zswitch = smask0(ji,jj) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice + IF( at_i(ji,jj) >= epsi10 ) THEN ; zswitch = smask0(ji,jj) + ELSE ; zswitch = 0._wp + ENDIF ! ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! zqld = smask0(ji,jj) * rDt_ice * & diff --git a/src/ICE/icetab.F90 b/src/ICE/icetab.F90 index 5db21b30578719112f6fe53c00bd9ab422e81c7e..db9575fcd28e30cf9868554b98461545fe73a740 100644 --- a/src/ICE/icetab.F90 +++ b/src/ICE/icetab.F90 @@ -55,11 +55,11 @@ CONTAINS ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls) ! - DO jn = 1, ndim1d - jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0 - jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0 - DO jl = 1, jpl - DO jk = 1, ipk + DO jl = 1, jpl + DO jk = 1, ipk + DO jn = 1, ndim1d + jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0 + jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0 tab1d(jn,jk,jl) = tab2d(jid,jjd,jk,jl) END DO END DO @@ -87,10 +87,10 @@ CONTAINS ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls) ! - DO jn = 1, ndim1d - jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0 - jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0 - DO jl = 1, jpl + DO jl = 1, jpl + DO jn = 1, ndim1d + jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0 + jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0 tab1d(jn,jl) = tab2d(jid,jjd,jl) END DO END DO @@ -145,11 +145,11 @@ CONTAINS ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls) ! - DO jn = 1, ndim1d - jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0 - jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0 - DO jl = 1, jpl - DO jk = 1, ipk + DO jl = 1, jpl + DO jk = 1, ipk + DO jn = 1, ndim1d + jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0 + jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0 tab2d(jid,jjd,jk,jl) = tab1d(jn,jk,jl) END DO END DO @@ -177,10 +177,10 @@ CONTAINS ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls) ! - DO jn = 1, ndim1d - jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0 - jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0 - DO jl = 1, jpl + DO jl = 1, jpl + DO jn = 1, ndim1d + jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0 + jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0 tab2d(jid,jjd,jl) = tab1d(jn,jl) END DO END DO diff --git a/src/ICE/icethd.F90 b/src/ICE/icethd.F90 index 592df9fa630ef8bae60cd21624ad9b9ed81d40a0..ddc4ddf2b6fdce2307bbf7f398c74e17f34c0b17 100644 --- a/src/ICE/icethd.F90 +++ b/src/ICE/icethd.F90 @@ -83,6 +83,7 @@ CONTAINS ! INTEGER :: ji, jj, jk, jl ! dummy loop indices !!------------------------------------------------------------------- + ! controls IF( ln_timing ) CALL timing_start('icethd') ! timing IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation @@ -201,20 +202,19 @@ CONTAINS !!------------------------------------------------------------------- INTEGER :: ji, jk ! dummy loop indices REAL(wp) :: ztmelts, zbbb, zccc ! local scalar - REAL(wp) :: zswitch !!------------------------------------------------------------------- ! Recover ice temperature DO jk = 1, nlay_i DO ji = 1, npti - ztmelts = -rTmlt * sz_i_1d(ji,jk) - ! Conversion q(S,T) -> T (second order equation) - zbbb = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus - zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) ) - t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi - - ! mask temperature - zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) - t_i_1d(ji,jk) = zswitch * t_i_1d(ji,jk) + ( 1._wp - zswitch ) * rt0 + IF( h_i_1d(ji) > 0._wp ) THEN + ztmelts = -rTmlt * sz_i_1d(ji,jk) + ! Conversion q(S,T) -> T (second order equation) + zbbb = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus + zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) ) + t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi + ELSE + t_i_1d(ji,jk) = rt0 + ENDIF END DO END DO ! @@ -232,7 +232,6 @@ CONTAINS REAL(wp) :: zhi_bef ! ice thickness before thermo REAL(wp) :: zdh_mel, zda_mel ! net melting REAL(wp) :: zvi, zvs ! ice/snow volumes - REAL(wp) :: zswitch !!----------------------------------------------------------------------- ! DO ji = 1, npti @@ -242,8 +241,7 @@ CONTAINS zvs = a_i_1d(ji) * h_s_1d(ji) ! lateral melting = concentration change zhi_bef = h_i_1d(ji) - zdh_mel - zswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) - zda_mel = zswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) + zda_mel = MAX( -a_i_1d(ji) , a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) ) a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel ) ! adjust thickness h_i_1d(ji) = zvi / a_i_1d(ji) diff --git a/src/ICE/icethd_dh.F90 b/src/ICE/icethd_dh.F90 index 0d7497f6284d96d09939ecafa953b9828d4ea835..5264b7566aef2a99676495afece415a797844e41 100644 --- a/src/ICE/icethd_dh.F90 +++ b/src/ICE/icethd_dh.F90 @@ -68,9 +68,6 @@ CONTAINS REAL(wp) :: ztmelts ! local scalar REAL(wp) :: zdum REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment - REAL(wp) :: zswi1 ! switch for computation of bottom salinity - REAL(wp) :: zswi12 ! switch for computation of bottom salinity - REAL(wp) :: zswi2 ! switch for computation of bottom salinity REAL(wp) :: zgrr ! bottom growth rate REAL(wp) :: zt_i_new ! bottom formation temperature REAL(wp) :: z1_rho ! 1/(rhos+rho0-rhoi) @@ -96,7 +93,7 @@ CONTAINS REAL(wp), DIMENSION(0:nlay_i+1) :: zh_i_old ! old thickness REAL(wp), DIMENSION(0:nlay_i+1) :: ze_i_old ! old enthalpy - REAL(wp) :: zswitch, zswitch_sal + REAL(wp) :: zswitch_sal INTEGER :: num_iter_max ! Heat conservation !!------------------------------------------------------------------ @@ -112,19 +109,22 @@ CONTAINS ! ! ============================================== ! IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN DO ji = 1, npti - zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice ) + zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice ) END DO ELSE DO ji = 1, npti - zdum = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji) - qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) - zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice ) + IF( t_su_1d(ji) >= rt0 ) THEN + qml_ice_1d(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji) + ELSE + qml_ice_1d(ji) = 0._wp + ENDIF + zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice ) END DO ENDIF ! DO ji = 1, npti - zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji) - zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice ) + zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji) + zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice ) END DO ! ! ========== ! ! ! Other init ! @@ -202,13 +202,12 @@ CONTAINS DO jk = 0, nlay_s IF( zh_s(jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN ! - zswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(jk) - epsi20 ) ) - zdum = - zswitch * zq_top(ji) / MAX( ze_s(jk), epsi20 ) ! thickness change - zdum = MAX( zdum , - zh_s(jk) ) ! bound melting - + zdum = - zq_top(ji) / MAX( ze_s(jk), epsi20 ) ! thickness change + zdum = MAX( zdum , - zh_s(jk) ) ! bound melting + hfx_snw_1d (ji) = hfx_snw_1d (ji) - ze_s(jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! heat used to melt snow(W.m-2, >0) wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! snow melting only = water into the ocean - + ! updates available heat + thickness dh_s_mlt(ji) = dh_s_mlt(ji) + zdum zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdum * ze_s(jk) ) @@ -324,8 +323,8 @@ CONTAINS ! record which layers have disappeared (for bottom melting) ! => icount=0 : no layer has vanished ! => icount=5 : 5 layers have vanished - zswitch = MAX( 0._wp , SIGN( 1._wp , - zh_i(jk) ) ) - icount(jk) = NINT( zswitch ) + IF( zh_i(jk) > 0._wp ) THEN ; icount(jk) = 0 + ELSE ; icount(jk) = 1 ; ENDIF END DO @@ -350,15 +349,12 @@ CONTAINS DO iter = 1, num_iter_max ! iterations ! New bottom ice salinity (Cox & Weeks, JGR88 ) - !--- zswi1 if dh/dt < 2.0e-8 - !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7 - !--- zswi2 if dh/dt > 3.6e-7 - zgrr = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_Dt_ice , epsi10 ) ) - zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) - zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) - zswi1 = 1. - zswi2 * zswi12 - zfracs = MIN( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & - & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 ) + zgrr = MIN( 1.0e-3_wp, MAX ( dh_i_bog(ji) * r1_Dt_ice , epsi10 ) ) + ! + IF ( zgrr < 2.0e-8_wp ) THEN ; zfracs = 0.12_wp + ELSEIF( zgrr >= 3.6e-7_wp ) THEN ; zfracs = MIN( 0.26_wp / ( 0.26_wp + 0.74_wp * EXP(-724300._wp*zgrr) ) , 0.5_wp ) + ELSE ; zfracs = MIN( 0.8925_wp + 0.0568_wp * LOG(100._wp*zgrr), 0.5_wp ) + ENDIF s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity @@ -605,7 +601,6 @@ CONTAINS ! INTEGER :: ji ! dummy loop indices INTEGER :: jk0, jk1 ! old/new layer indices - REAL(wp) :: zswitch ! REAL(wp), DIMENSION(0:nlay_s+1) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces REAL(wp), DIMENSION(0:nlay_s) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces @@ -650,8 +645,7 @@ CONTAINS ! new enthalpies DO jk1 = 1, nlay_s - zswitch = MAX( 0._wp , SIGN( 1._wp , zhnew - epsi20 ) ) - snw_ent(jk1) = zswitch * ( zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 ) + snw_ent(jk1) = MAX( 0._wp, zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 ) ! max for roundoff error END DO @@ -687,7 +681,6 @@ CONTAINS ! INTEGER :: ji ! dummy loop indices INTEGER :: jk0, jk1 ! old/new layer indices - REAL(wp) :: zswitch ! REAL(wp), DIMENSION(0:nlay_i+2) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces REAL(wp), DIMENSION(0:nlay_i) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces @@ -732,8 +725,7 @@ CONTAINS ! new enthalpies DO jk1 = 1, nlay_i - zswitch = MAX( 0._wp , SIGN( 1._wp , zhnew - epsi20 ) ) - ice_ent1(jk1) = zswitch * MAX( 0._wp, zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 ) ! max for roundoff error + ice_ent1(jk1) = MAX( 0._wp, zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 ) ! max for roundoff error END DO ! --- diag error on heat remapping --- ! diff --git a/src/ICE/icethd_do.F90 b/src/ICE/icethd_do.F90 index 09cf8f65a8305631e6099840bcef61c7057a1a02..90e0cee2bfb606cab3ddbe26d7184071ccee394c 100644 --- a/src/ICE/icethd_do.F90 +++ b/src/ICE/icethd_do.F90 @@ -78,7 +78,6 @@ CONTAINS REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean) ! INTEGER :: jcat ! indexes of categories where new ice grows - REAL(wp) :: zswinew ! switch for new ice or not ! REAL(wp) :: zv_newfra REAL(wp) :: zv_newice ! volume of accreted ice @@ -88,7 +87,6 @@ CONTAINS REAL(wp) :: zdv_res ! residual volume in case of excessive heat budget REAL(wp) :: zda_res ! residual area in case of excessive heat budget REAL(wp) :: zv_frazb ! accretion of frazil ice at the ice bottom - REAL(wp) :: zswitch ! REAL(wp), DIMENSION(jpl) :: zv_b ! old volume of ice in category jl REAL(wp), DIMENSION(jpl) :: za_b ! old area of ice in category jl @@ -206,10 +204,12 @@ CONTAINS sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice * rhoi * zs_newice(ji) * r1_Dt_ice ! A fraction fraz_frac of frazil ice is accreted at the ice bottom - zswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) ) - zv_frazb = zfraz_frac_1d(ji) * zswitch * zv_newice - zv_newice = ( 1._wp - zfraz_frac_1d(ji) * zswitch ) * zv_newice - + IF( at_i_1d(ji) > 0._wp ) THEN + zv_frazb = zfraz_frac_1d(ji) * zv_newice + zv_newice = ( 1._wp - zfraz_frac_1d(ji) ) * zv_newice + ELSE + zv_frazb = 0._wp + ENDIF ! --- Area of new ice --- ! za_newice = zv_newice / zh_newice(ji) @@ -242,15 +242,11 @@ CONTAINS ! Heat content jl = jcat ! categroy in which new ice is put - zswinew = MAX( 0._wp , SIGN( 1._wp , - za_b(jl) ) ) ! 0 if old ice - - DO jk = 1, nlay_i - jl = jcat - zswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) ) - e_i_2d(ji,jk,jl) = zswinew * ze_newice + & - & ( 1.0 - zswinew ) * ( ze_newice * zv_newice + e_i_2d(ji,jk,jl) * zv_b(jl) ) & - & * zswitch / MAX( v_i_2d(ji,jl), epsi20 ) - END DO + IF( za_b(jl) > 0._wp ) THEN + e_i_2d(ji,:,jl) = ( ze_newice * zv_newice + e_i_2d(ji,:,jl) * zv_b(jl) ) / MAX( v_i_2d(ji,jl), epsi20 ) + ELSE + e_i_2d(ji,:,jl) = ze_newice + ENDIF ! --- bottom ice growth + ice enthalpy remapping --- ! DO jl = 1, jpl @@ -264,9 +260,12 @@ CONTAINS END DO ! new volumes including lateral/bottom accretion + residual - zswitch = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) ) - zv_newfra = zswitch * ( zdv_res + zv_frazb ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 ) - a_i_2d(ji,jl) = zswitch * a_i_2d(ji,jl) + IF( at_i_1d(ji) >= epsi20 ) THEN + zv_newfra = ( zdv_res + zv_frazb ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 ) + ELSE + zv_newfra = 0._wp + a_i_2d(ji,jl) = 0._wp + ENDIF v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra ! for remapping zh_i_old(nlay_i+1) = zv_newfra @@ -323,7 +322,6 @@ CONTAINS INTEGER :: ji, jj ! dummy loop indices INTEGER :: iter REAL(wp) :: zvfrx, zvgx, ztaux, zf, ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp - REAL(wp) :: zswitch REAL(wp), PARAMETER :: zcai = 1.4e-3_wp ! ice-air drag (clem: should be dependent on coupling/forcing used) REAL(wp), PARAMETER :: zhicrit = 0.04_wp ! frazil ice thickness REAL(wp), PARAMETER :: zsqcd = 1.0_wp / SQRT( 1.3_wp * zcai ) ! 1/SQRT(airdensity*drag) @@ -359,20 +357,25 @@ CONTAINS ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) ! -- Frazil ice velocity -- ! - zswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) - zvfrx = zswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) - zvfry = zswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) - + IF( ztenagm >= epsi10 ) THEN + zvfrx = zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) + zvfry = zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) + ELSE + zvfrx = 0._wp + zvfry = 0._wp + ENDIF ! -- Pack ice velocity -- ! - zvgx = ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp - zvgy = ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp - - ! -- Relative frazil/pack ice velocity -- ! - zswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) - zvrel2 = MAX( (zvfrx - zvgx)*(zvfrx - zvgx) + (zvfry - zvgy)*(zvfry - zvgy), 0.15_wp*0.15_wp ) * zswitch - - ! -- fraction of frazil ice -- ! - fraz_frac(ji,jj) = zswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz + zvgx = ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp + zvgy = ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp + + ! -- Relative frazil/pack ice velocity & fraction of frazil ice-- ! + IF( at_i(ji,jj) >= epsi10 ) THEN + zvrel2 = MAX( (zvfrx - zvgx)*(zvfrx - zvgx) + (zvfry - zvgy)*(zvfry - zvgy), 0.15_wp*0.15_wp ) + fraz_frac(ji,jj) = ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz + ELSE + zvrel2 = 0._wp + fraz_frac(ji,jj) = 0._wp + ENDIF ! -- new ice thickness (iterative loop) -- ! ht_i_new(ji,jj) = zhicrit + ( zhicrit + 0.1_wp ) & @@ -429,7 +432,6 @@ CONTAINS ! INTEGER :: ji ! dummy loop indices INTEGER :: jk0, jk1 ! old/new layer indices - REAL(wp) :: zswitch ! REAL(wp), DIMENSION(0:nlay_i+2) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces REAL(wp), DIMENSION(0:nlay_i) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces @@ -474,8 +476,7 @@ CONTAINS ! new enthalpies DO jk1 = 1, nlay_i - zswitch = MAX( 0._wp , SIGN( 1._wp , zhnew - epsi20 ) ) - ice_ent2(jk1) = zswitch * MAX( 0._wp, zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 ) ! max for roundoff error + ice_ent2(jk1) = MAX( 0._wp, zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 ) ! max for roundoff error END DO ! --- diag error on heat remapping --- ! diff --git a/src/ICE/icethd_zdf_bl99.F90 b/src/ICE/icethd_zdf_bl99.F90 index 979aebc30b189ef3caa310e563de7eedb0d8c521..c68662595d3c33c172e725245536cde4038a598c 100644 --- a/src/ICE/icethd_zdf_bl99.F90 +++ b/src/ICE/icethd_zdf_bl99.F90 @@ -206,18 +206,19 @@ CONTAINS ! 2) Radiation !------------- ! --- Transmission/absorption of solar radiation in the ice --- ! - DO ji = 1, npti - ! - zradtr_s(ji,0) = qtr_ice_top_1d(ji) - DO jk = 1, nlay_s + zradtr_s(1:npti,0) = qtr_ice_top_1d(1:npti) + DO jk = 1, nlay_s + DO ji = 1, npti ! ! radiation transmitted below the layer-th snow layer zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s(ji) * MAX( 0._wp, zh_s(ji) * REAL(jk) - zhs_ssl ) ) ! ! radiation absorbed by the layer-th snow layer zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) END DO - ! - zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * za_s_fra(ji) + qtr_ice_top_1d(ji) * ( 1._wp - za_s_fra(ji) ) - DO jk = 1, nlay_i + END DO + ! + zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * za_s_fra(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - za_s_fra(1:npti) ) + DO jk = 1, nlay_i + DO ji = 1, npti ! ! radiation transmitted below the layer-th ice layer zradtr_i(ji,jk) = za_s_fra(ji) * zradtr_s(ji,nlay_s) & ! part covered by snow & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zh_min ) ) & @@ -226,9 +227,10 @@ CONTAINS ! ! radiation absorbed by the layer-th ice layer zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) END DO + END DO ! + DO ji = 1, npti qtr_ice_bot_1d(ji) = zradtr_i(ji,nlay_i) ! record radiation transmitted below the ice - ! END DO ! ! @@ -252,7 +254,9 @@ CONTAINS DO ji = 1, npti ztcond_i_cp(ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) - DO jk = 1, nlay_i-1 + END DO + DO jk = 1, nlay_i-1 + DO ji = 1, npti ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & & MIN( -epsi10, 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) END DO @@ -265,7 +269,9 @@ CONTAINS & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) - DO jk = 1, nlay_i-1 + END DO + DO jk = 1, nlay_i-1 + DO ji = 1, npti ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & & MIN( -epsi10, 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) & & - 0.011_wp * ( 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) @@ -299,10 +305,11 @@ CONTAINS ! ENDIF ! + !----------------- + ! 4) kappa factors + !----------------- DO ji = 1, npti - !----------------- - ! 4) kappa factors - !----------------- + ! IF ( .NOT. l_T_converged(ji) ) THEN ! !--- Snow @@ -325,20 +332,28 @@ CONTAINS ! If there is snow then use the same snow-ice interface conductivity for the top layer of ice IF( h_s_1d(ji) > 0._wp ) zkappa_i(ji,0) = zkappa_s(ji,nlay_s) ! Snow-ice interface ! - !-------------------------------------- - ! 5) Sea ice specific heat, eta factors - !-------------------------------------- - DO jk = 1, nlay_i - zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) - zeta_i(ji,jk) = rDt_ice * r1_rhoi * z1_h_i(ji) / zcpi - END DO - DO jk = 1, nlay_s - zeta_s(ji,jk) = rDt_ice * r1_rhos * r1_rcpi * z1_h_s(ji) - END DO - ! ENDIF ! END DO + + !-------------------------------------- + ! 5) Sea ice specific heat, eta factors + !-------------------------------------- + DO jk = 1, nlay_i + DO ji = 1, npti + IF ( .NOT. l_T_converged(ji) ) THEN + zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) + zeta_i(ji,jk) = rDt_ice * r1_rhoi * z1_h_i(ji) / zcpi + ENDIF + END DO + END DO + ! + DO jk = 1, nlay_s + DO ji = 1, npti + IF ( .NOT. l_T_converged(ji) ) & + & zeta_s(ji,jk) = rDt_ice * r1_rhos * r1_rcpi * z1_h_s(ji) + END DO + END DO ! ! IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN @@ -365,14 +380,9 @@ CONTAINS ! zfnet = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar ! - ! ! before temperatures - DO jk = 1, nlay_i - ztib(jk) = t_i_1d(ji,jk) - ENDDO - DO jk = 1, nlay_s - ztsb(jk) = t_s_1d(ji,jk) - ENDDO + ztib(:) = t_i_1d(ji,:) + ztsb(:) = t_s_1d(ji,:) ! !---------------------------- ! 7) tridiagonal system terms @@ -612,12 +622,8 @@ CONTAINS ! IF ( .NOT. l_T_converged(ji) ) THEN ! before temperatures - DO jk = 1, nlay_i - ztib(jk) = t_i_1d(ji,jk) - ENDDO - DO jk = 1, nlay_s - ztsb(jk) = t_s_1d(ji,jk) - ENDDO + ztib(:) = t_i_1d(ji,:) + ztsb(:) = t_s_1d(ji,:) ! !---------------------------- ! 7) tridiagonal system terms diff --git a/src/ICE/iceupdate.F90 b/src/ICE/iceupdate.F90 index dab6f7a62ecef608bd841dc27672cdaecde11418..f9806bd58e074c4aa43555f4d129c63f0965c0e0 100644 --- a/src/ICE/iceupdate.F90 +++ b/src/ICE/iceupdate.F90 @@ -177,7 +177,7 @@ CONTAINS wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj) ! total mass flux at the ocean/ice interface - fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model + fwfice(ji,jj) = wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_pnd(ji,jj) + wfx_err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux ! Salt flux at the ocean surface diff --git a/src/ICE/icevar.F90 b/src/ICE/icevar.F90 index fe8a20999f2669a5b6792db115a95f7f4928bd7b..03d1fd6d28e5c7f8f1c1a2047f03237a0f56ee1f 100644 --- a/src/ICE/icevar.F90 +++ b/src/ICE/icevar.F90 @@ -114,8 +114,7 @@ CONTAINS ! ! >1 state variables + others ! INTEGER :: ji, jj, jk, jl ! dummy loop indices - REAL(wp) :: z1_vt_i, z1_vt_s - REAL(wp), DIMENSION(A2D(0)) :: z1_at_i + REAL(wp) :: z1_vt_i, z1_vt_s, z1_at_i !!------------------------------------------------------------------- ! ! full arrays: vt_i, vt_s, at_i, vt_ip, vt_il, at_ip @@ -141,11 +140,9 @@ CONTAINS ! !!GS: tm_su always needed by ABL over sea-ice IF( at_i(ji,jj) <= epsi20 ) THEN - z1_at_i(ji,jj) = 0._wp tm_su (ji,jj) = rt0 ELSE - z1_at_i(ji,jj) = 1._wp / at_i(ji,jj) - tm_su (ji,jj) = SUM( t_su(ji,jj,:) * a_i(ji,jj,:) ) * z1_at_i(ji,jj) + tm_su (ji,jj) = SUM( t_su(ji,jj,:) * a_i(ji,jj,:) ) / at_i(ji,jj) ENDIF END_2D ! @@ -154,34 +151,76 @@ CONTAINS IF( kn > 1 ) THEN ! DO_2D( 0, 0, 0, 0 ) + IF( at_i(ji,jj) > epsi20 ) THEN ; z1_at_i = 1._wp / at_i(ji,jj) + ELSE ; z1_at_i = 0._wp + ENDIF IF( vt_i(ji,jj) > epsi20 ) THEN ; z1_vt_i = 1._wp / vt_i(ji,jj) ELSE ; z1_vt_i = 0._wp ENDIF - IF( vt_s(ji,jj) > epsi20 ) THEN ; z1_vt_s = 1._wp / vt_s(ji,jj) - ELSE ; z1_vt_s = 0._wp - ENDIF ! mean ice/snow thickness - hm_i(ji,jj) = vt_i(ji,jj) * z1_at_i(ji,jj) - hm_s(ji,jj) = vt_s(ji,jj) * z1_at_i(ji,jj) + hm_i(ji,jj) = vt_i(ji,jj) * z1_at_i + hm_s(ji,jj) = vt_s(ji,jj) * z1_at_i ! ! mean temperature (K), salinity and age - tm_si(ji,jj) = SUM( t_si(ji,jj,:) * a_i(ji,jj,:) ) * z1_at_i(ji,jj) - om_i (ji,jj) = SUM( oa_i(ji,jj,:) ) * z1_at_i(ji,jj) + tm_si(ji,jj) = SUM( t_si(ji,jj,:) * a_i(ji,jj,:) ) * z1_at_i + om_i (ji,jj) = SUM( oa_i(ji,jj,:) ) * z1_at_i sm_i (ji,jj) = st_i(ji,jj) * z1_vt_i + END_2D ! - tm_i(ji,jj) = 0._wp - tm_s(ji,jj) = 0._wp - DO jl = 1, jpl - DO jk = 1, nlay_i - tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * t_i (ji,jj,jk,jl) * v_i(ji,jj,jl) * z1_vt_i - END DO - DO jk = 1, nlay_s - tm_s(ji,jj) = tm_s(ji,jj) + r1_nlay_s * t_s (ji,jj,jk,jl) * v_s(ji,jj,jl) * z1_vt_s - END DO - END DO + tm_i(:,:) = 0._wp + tm_s(:,:) = 0._wp + DO jl = 1, jpl + DO_3D( 0, 0, 0, 0, 1, nlay_i ) + IF( vt_i(ji,jj) > epsi20 ) THEN + tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * t_i (ji,jj,jk,jl) * v_i(ji,jj,jl) / vt_i(ji,jj) + ELSE + tm_i(ji,jj) = rt0 + ENDIF + END_3D + END DO + DO jl = 1, jpl + DO_3D( 0, 0, 0, 0, 1, nlay_s ) + IF( vt_s(ji,jj) > epsi20 ) THEN + tm_s(ji,jj) = tm_s(ji,jj) + r1_nlay_s * t_s (ji,jj,jk,jl) * v_s(ji,jj,jl) / vt_s(ji,jj) + ELSE + tm_s(ji,jj) = rt0 + ENDIF + END_3D + END DO ! - END_2D +!!$ DO_2D( 0, 0, 0, 0 ) +!!$ IF( at_i(ji,jj) > epsi20 ) THEN ; z1_at_i = 1._wp / at_i(ji,jj) +!!$ ELSE ; z1_at_i = 0._wp +!!$ ENDIF +!!$ IF( vt_i(ji,jj) > epsi20 ) THEN ; z1_vt_i = 1._wp / vt_i(ji,jj) +!!$ ELSE ; z1_vt_i = 0._wp +!!$ ENDIF +!!$ IF( vt_s(ji,jj) > epsi20 ) THEN ; z1_vt_s = 1._wp / vt_s(ji,jj) +!!$ ELSE ; z1_vt_s = 0._wp +!!$ ENDIF +!!$ +!!$ ! mean ice/snow thickness +!!$ hm_i(ji,jj) = vt_i(ji,jj) * z1_at_i +!!$ hm_s(ji,jj) = vt_s(ji,jj) * z1_at_i +!!$ ! +!!$ ! mean temperature (K), salinity and age +!!$ tm_si(ji,jj) = SUM( t_si(ji,jj,:) * a_i(ji,jj,:) ) * z1_at_i +!!$ om_i (ji,jj) = SUM( oa_i(ji,jj,:) ) * z1_at_i +!!$ sm_i (ji,jj) = st_i(ji,jj) * z1_vt_i +!!$ ! +!!$ tm_i(ji,jj) = 0._wp +!!$ tm_s(ji,jj) = 0._wp +!!$ DO jl = 1, jpl +!!$ DO jk = 1, nlay_i +!!$ tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * t_i (ji,jj,jk,jl) * v_i(ji,jj,jl) * z1_vt_i +!!$ END DO +!!$ DO jk = 1, nlay_s +!!$ tm_s(ji,jj) = tm_s(ji,jj) + r1_nlay_s * t_s (ji,jj,jk,jl) * v_s(ji,jj,jl) * z1_vt_s +!!$ END DO +!!$ END DO +!!$ ! +!!$ END_2D ! put rt0 where there is no ice WHERE( at_i(A2D(0)) <= epsi20 ) tm_si(:,:) = rt0 @@ -215,12 +254,11 @@ CONTAINS INTEGER :: ji, jj, jk, jl ! dummy loop indices REAL(wp) :: ze_i ! local scalars REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - - - REAL(wp) :: zhmax, z1_zhmax ! - - + REAL(wp) :: zhmax, z1_hmax ! - - REAL(wp) :: zlay_i, zlay_s ! - - REAL(wp), PARAMETER :: zhl_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation REAL(wp), PARAMETER :: zhl_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation - REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i - REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z1_a_ip + REAL(wp) :: z1_hl, z1_a_i, z1_a_ip REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: za_s_fra !!------------------------------------------------------------------- @@ -232,97 +270,109 @@ CONTAINS !--------------------------------------------------------------- ! Ice thickness, snow thickness, ice salinity, ice age and ponds !--------------------------------------------------------------- - ! !--- inverse of the ice area - WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) - ELSEWHERE ; z1_a_i(:,:,:) = 0._wp - END WHERE - ! - WHERE( v_i(:,:,:) > epsi20 ) ; z1_v_i(:,:,:) = 1._wp / v_i(:,:,:) - ELSEWHERE ; z1_v_i(:,:,:) = 0._wp - END WHERE ! - ! !--- ice thickness - h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) + ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area + ! clem: if a>1 then do something + zhmax = hi_max(jpl) + z1_hmax = 1._wp / hi_max(jpl) + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) + IF( v_i(ji,jj,jpl) > ( zhmax*a_i(ji,jj,jpl) ) ) a_i(ji,jj,jpl) = MIN( 1._wp, v_i(ji,jj,jpl) * z1_hmax ) + END_2D + + DO jl = 1, jpl + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) + ! !--- inverse of the ice area + IF( a_i(ji,jj,jl) > epsi20 ) THEN ; z1_a_i = 1._wp / a_i(ji,jj,jl) + ELSE ; z1_a_i = 0._wp + ENDIF + ! !--- ice thickness + h_i(ji,jj,jl) = v_i (ji,jj,jl) * z1_a_i + ! !--- snow thickness + h_s(ji,jj,jl) = v_s (ji,jj,jl) * z1_a_i + ! !--- ice age + o_i(ji,jj,jl) = oa_i(ji,jj,jl) * z1_a_i + ! + END_2D + ENDDO + ! !--- salinity (with a minimum value imposed everywhere) + IF( nn_icesal == 2 ) THEN + WHERE( v_i(:,:,:) > epsi20 ) ; s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) / v_i(:,:,:) ) ) + ELSEWHERE ; s_i(:,:,:) = rn_simin + END WHERE + ENDIF + CALL ice_var_salprof ! salinity profile - zhmax = hi_max(jpl) - z1_zhmax = 1._wp / hi_max(jpl) - WHERE( h_i(:,:,jpl) > zhmax ) ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area - h_i (:,:,jpl) = zhmax - a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax - z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) - END WHERE - ! !--- snow thickness - h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:) - ! !--- ice age - o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) - ! !--- pond and lid thickness IF( kn == 1 .OR. ln_pnd ) THEN - ALLOCATE( z1_a_ip(jpi,jpj,jpl), za_s_fra(A2D(0),jpl) ) - ! - WHERE( a_ip(:,:,:) > epsi20 ) ; z1_a_ip(:,:,:) = 1._wp / a_ip(:,:,:) - ELSEWHERE ; z1_a_ip(:,:,:) = 0._wp - END WHERE + ALLOCATE( za_s_fra(A2D(0),jpl) ) ! - h_ip(:,:,:) = v_ip(:,:,:) * z1_a_ip(:,:,:) - h_il(:,:,:) = v_il(:,:,:) * z1_a_ip(:,:,:) - ! !--- melt pond effective area (used for albedo) - a_ip_frac(:,:,:) = a_ip(A2D(0),:) * z1_a_i(A2D(0),:) - WHERE ( h_il(A2D(0),:) <= zhl_min ) ; a_ip_eff(:,:,:) = a_ip_frac(:,:,:) ! lid is very thin. Expose all the pond - ELSEWHERE( h_il(A2D(0),:) >= zhl_max ) ; a_ip_eff(:,:,:) = 0._wp ! lid is very thick. Cover all the pond up with ice and snow - ELSEWHERE ; a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * & ! lid is in between. Expose part of the pond - & ( zhl_max - h_il(A2D(0),:) ) / ( zhl_max - zhl_min ) - END WHERE + z1_hl = 1._wp / ( zhl_max - zhl_min ) + DO jl = 1, jpl + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) + IF( a_ip(ji,jj,jl) > epsi20 ) THEN ; z1_a_ip = 1._wp / a_ip(ji,jj,jl) + ELSE ; z1_a_ip = 0._wp + ENDIF + ! !--- pond and lid thickness + h_ip(ji,jj,jl) = v_ip(ji,jj,jl) * z1_a_ip + h_il(ji,jj,jl) = v_il(ji,jj,jl) * z1_a_ip + END_2D + ! !--- melt pond effective area (used for albedo) + DO_2D( 0, 0, 0, 0 ) + IF( a_i(ji,jj,jl) > epsi20 ) THEN ; a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) + ELSE ; a_ip_frac(ji,jj,jl) = 0._wp + ENDIF + IF ( h_il(ji,jj,jl) <= zhl_min ) THEN ; a_ip_eff(ji,jj,jl) = a_ip_frac(ji,jj,jl) ! lid is very thin. Expose all the pond + ELSEIF( h_il(ji,jj,jl) >= zhl_max ) THEN ; a_ip_eff(ji,jj,jl) = 0._wp ! lid is very thick. Cover all the pond up with ice and snow + ELSE ; a_ip_eff(ji,jj,jl) = a_ip_frac(ji,jj,jl) * & ! lid is in between. Expose part of the pond + & ( zhl_max - h_il(ji,jj,jl) ) * z1_hl + ENDIF + ! + END_2D + ENDDO ! - CALL ice_var_snwfra( h_s(A2D(0),:), za_s_fra(:,:,:) ) ! calculate ice fraction covered by snow + CALL ice_var_snwfra( h_s(A2D(0),:), za_s_fra(:,:,:) ) ! calculate ice fraction covered by snow a_ip_eff(:,:,:) = MIN( a_ip_eff(:,:,:), 1._wp - za_s_fra(:,:,:) ) ! make sure (a_ip_eff + a_s_fra) <= 1 ! - DEALLOCATE( z1_a_ip, za_s_fra ) - ENDIF - ! - ! !--- salinity (with a minimum value imposed everywhere) - IF( nn_icesal == 2 ) THEN - WHERE( v_i(:,:,:) > epsi20 ) ; s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) ) - ELSEWHERE ; s_i(:,:,:) = rn_simin - END WHERE + DEALLOCATE( za_s_fra ) ENDIF - CALL ice_var_salprof ! salinity profile - + !------------------- ! Ice temperature [K] (with a minimum value (rt0 - 100.)) !------------------- - zlay_i = REAL( nlay_i , wp ) ! number of layers + zlay_i = REAL( nlay_i , wp ) ! number of layers DO jl = 1, jpl - DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) - IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area - DO jk = 1, nlay_i + DO jk = 1, nlay_i + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) + IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area ! - ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] - ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C] + ze_i = e_i (ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] + ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C] ! Conversion q(S,T) -> T (second order equation) zbbb = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) ) t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts ! - END DO - ELSE !--- no ice - DO jk = 1, nlay_i + ELSE !--- no ice t_i(ji,jj,jk,jl) = rt0 - END DO - ENDIF - END_2D + ENDIF + END_2D + END DO END DO !-------------------- ! Snow temperature [K] (with a minimum value (rt0 - 100.)) !-------------------- zlay_s = REAL( nlay_s , wp ) - DO jk = 1, nlay_s - WHERE( v_s(:,:,:) > epsi20 ) !--- icy area - t_s(:,:,jk,:) = rt0 + MAX( -100._wp , & - & MIN( r1_rcpi * ( -r1_rhos * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + rLfus ) , 0._wp ) ) - ELSEWHERE !--- no ice - t_s(:,:,jk,:) = rt0 - END WHERE + DO jl = 1, jpl + DO jk = 1, nlay_s + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) + IF ( v_s(ji,jj,jl) > epsi20 ) THEN !--- icy area + t_s(ji,jj,jk,jl) = rt0 + MAX( -100._wp , & + & MIN( r1_rcpi*( -r1_rhos*( e_s(ji,jj,jk,jl) / v_s(ji,jj,jl) * zlay_s ) + rLfus ) , 0._wp ) ) + ELSE !--- no ice + t_s(ji,jj,jk,jl) = rt0 + ENDIF + END_2D + END DO END DO ! ! integrated values @@ -394,26 +444,26 @@ CONTAINS ! !---------------------------------------------! z1_dS = 1._wp / ( zsi1 - zsi0 ) ! - DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) - DO jl = 1, jpl - ! ! Slope of the linear profile - IF( h_i(ji,jj,jl) > epsi20 ) THEN ; z_slope_s = 2._wp * s_i(ji,jj,jl) / h_i(ji,jj,jl) - ELSE ; z_slope_s = 0._wp - ENDIF - ! - zalpha = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) ) - ! ! force a constant profile when SSS too low (Baltic Sea) - IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha = 0._wp - ! - ! Computation of the profile - DO jk = 1, nlay_i + DO jl = 1, jpl + DO jk = 1, nlay_i + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) + ! ! Slope of the linear profile + IF( h_i(ji,jj,jl) > epsi20 ) THEN ; z_slope_s = 2._wp * s_i(ji,jj,jl) / h_i(ji,jj,jl) + ELSE ; z_slope_s = 0._wp + ENDIF + ! + zalpha = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) ) + ! ! force a constant profile when SSS too low (Baltic Sea) + IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha = 0._wp + ! + ! Computation of the profile ! ! linear profile with 0 surface value zs0 = z_slope_s * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i zs = zalpha * zs0 + ( 1._wp - zalpha ) * s_i(ji,jj,jl) ! weighting the profile sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) ) - ENDDO + END_2D ENDDO - END_2D + ENDDO ! ! !-------------------------------------------! CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile @@ -469,19 +519,19 @@ CONTAINS ! !---------------------------------------------! z1_dS = 1._wp / ( zsi1 - zsi0 ) ! - DO ji = 1, npti - ! ! Slope of the linear profile - IF( h_i_1d(ji) > epsi20 ) THEN ; z_slope_s = 2._wp * s_i_1d(ji) / h_i_1d(ji) - ELSE ; z_slope_s = 0._wp - ENDIF - ! - zalpha = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) ) - ! ! force a constant profile when SSS too low (Baltic Sea) - IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha = 0._wp - ! - ! - ! Computation of the profile - DO jk = 1, nlay_i + DO jk = 1, nlay_i + DO ji = 1, npti + ! ! Slope of the linear profile + IF( h_i_1d(ji) > epsi20 ) THEN ; z_slope_s = 2._wp * s_i_1d(ji) / h_i_1d(ji) + ELSE ; z_slope_s = 0._wp + ENDIF + ! + zalpha = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) ) + ! ! force a constant profile when SSS too low (Baltic Sea) + IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha = 0._wp + ! + ! + ! Computation of the profile ! ! linear profile with 0 surface value zs0 = z_slope_s * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i zs = zalpha * zs0 + ( 1._wp - zalpha ) * s_i_1d(ji) @@ -516,66 +566,78 @@ CONTAINS !! ** Purpose : Remove too small sea ice areas and correct fluxes !!------------------------------------------------------------------- INTEGER :: ji, jj, jl, jk ! dummy loop indices - REAL(wp), DIMENSION(A2D(0)) :: zswitch + REAL(wp) :: zsmall !!------------------------------------------------------------------- ! - DO jl = 1, jpl !== loop over the categories ==! - ! - WHERE( a_i(A2D(0),jl) > epsi10 ) ; h_i(A2D(0),jl) = v_i(A2D(0),jl) / a_i(A2D(0),jl) - ELSEWHERE ; h_i(A2D(0),jl) = 0._wp - END WHERE - ! - WHERE( a_i(A2D(0),jl) < epsi10 .OR. v_i(A2D(0),jl) < epsi10 .OR. h_i(A2D(0),jl) < epsi10 ) ; zswitch(:,:) = 0._wp - ELSEWHERE ; zswitch(:,:) = 1._wp - END WHERE - ! - DO_2D( 0, 0, 0, 0 ) - !----------------------------------------------------------------- - ! Zap ice energy and use ocean heat to melt ice - !----------------------------------------------------------------- - DO jk = 1, nlay_i - ! update exchanges with ocean - hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 - e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj) - t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) - ENDDO + WHERE( a_i(A2D(0),:) > epsi10 ) ; h_i(A2D(0),:) = v_i(A2D(0),:) / a_i(A2D(0),:) + ELSEWHERE ; h_i(A2D(0),:) = 0._wp + END WHERE + ! + !----------------------------------------------------------------- + ! Zap ice energy and use ocean heat to melt ice + !----------------------------------------------------------------- + DO jl = 1, jpl + DO_3D( 0, 0, 0, 0, 1, nlay_i ) ! - DO jk = 1, nlay_s - ! update exchanges with ocean - hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 - e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj) - t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) - ENDDO + zsmall = MIN( a_i(ji,jj,jl), v_i(ji,jj,jl), h_i(ji,jj,jl) ) ! - !----------------------------------------------------------------- - ! zap ice and snow volume, add water and salt to ocean - !----------------------------------------------------------------- - ! update exchanges with ocean - sfx_res(ji,jj) = sfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_Dt_ice - wfx_res(ji,jj) = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_Dt_ice - wfx_res(ji,jj) = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_Dt_ice - wfx_res(ji,jj) = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * ( v_ip(ji,jj,jl)+v_il(ji,jj,jl) ) * rhow * r1_Dt_ice + IF( zsmall < epsi10 ) THEN + ! update exchanges with ocean + hfx_res(ji,jj) = hfx_res(ji,jj) - e_i(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 + e_i(ji,jj,jk,jl) = 0._wp + t_i(ji,jj,jk,jl) = rt0 + ENDIF + END_3D + ENDDO + + DO jl = 1, jpl + DO_3D( 0, 0, 0, 0, 1, nlay_s ) ! - a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) - v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj) - v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj) - t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + ( sst_m(ji,jj) + rt0 ) * ( 1._wp - zswitch(ji,jj) ) - oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj) - sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj) + zsmall = MIN( a_i(ji,jj,jl), v_i(ji,jj,jl), h_i(ji,jj,jl) ) ! - h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj) - h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj) + IF( zsmall < epsi10 ) THEN + ! update exchanges with ocean + hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 + e_s(ji,jj,jk,jl) = 0._wp + t_s(ji,jj,jk,jl) = rt0 + ENDIF + END_3D + ENDDO + ! + !----------------------------------------------------------------- + ! zap ice and snow volume, add water and salt to ocean + !----------------------------------------------------------------- + DO jl = 1, jpl + DO_2D( 0, 0, 0, 0 ) ! - a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) - v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) - v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj) - h_ip (ji,jj,jl) = h_ip (ji,jj,jl) * zswitch(ji,jj) - h_il (ji,jj,jl) = h_il (ji,jj,jl) * zswitch(ji,jj) + zsmall = MIN( a_i(ji,jj,jl), v_i(ji,jj,jl), h_i(ji,jj,jl) ) ! + IF( zsmall < epsi10 ) THEN + ! update exchanges with ocean + sfx_res(ji,jj) = sfx_res(ji,jj) + sv_i(ji,jj,jl) * rhoi * r1_Dt_ice + wfx_res(ji,jj) = wfx_res(ji,jj) + v_i (ji,jj,jl) * rhoi * r1_Dt_ice + wfx_res(ji,jj) = wfx_res(ji,jj) + v_s (ji,jj,jl) * rhos * r1_Dt_ice + wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_ip(ji,jj,jl)+v_il(ji,jj,jl) ) * rhow * r1_Dt_ice + ! + a_i (ji,jj,jl) = 0._wp + v_i (ji,jj,jl) = 0._wp + v_s (ji,jj,jl) = 0._wp + t_su (ji,jj,jl) = sst_m(ji,jj) + rt0 + oa_i (ji,jj,jl) = 0._wp + sv_i (ji,jj,jl) = 0._wp + ! + h_i (ji,jj,jl) = 0._wp + h_s (ji,jj,jl) = 0._wp + ! + a_ip (ji,jj,jl) = 0._wp + v_ip (ji,jj,jl) = 0._wp + v_il (ji,jj,jl) = 0._wp + h_ip (ji,jj,jl) = 0._wp + h_il (ji,jj,jl) = 0._wp + ENDIF END_2D - ! END DO - + ! to be sure that at_i is the sum of a_i(jl) at_i (A2D(0)) = SUM( a_i (A2D(0),:), dim=3 ) vt_i (A2D(0)) = SUM( v_i (A2D(0),:), dim=3 ) @@ -625,32 +687,35 @@ CONTAINS z1_dt = 1._wp / pdt ! - DO jl = 1, jpl !== loop over the categories ==! - ! - ! make sure a_i=0 where v_i<=0 - WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp - + ! make sure a_i=0 where v_i<=0 + WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp + + !---------------------------------------- + ! zap ice energy and send it to the ocean + !---------------------------------------- + DO jl = 1, jpl + DO_3D( ihls, ihls, ihls, ihls, 1, nlay_i ) + IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN + zhfx_res(ji,jj) = zhfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 + pe_i(ji,jj,jk,jl) = 0._wp + ENDIF + END_3D + ENDDO + ! + DO jl = 1, jpl + DO_3D( ihls, ihls, ihls, ihls, 1, nlay_s ) + IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN + zhfx_res(ji,jj) = zhfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 + pe_s(ji,jj,jk,jl) = 0._wp + ENDIF + END_3D + ENDDO + ! + !----------------------------------------------------- + ! zap ice and snow volume, add water and salt to ocean + !----------------------------------------------------- + DO jl = 1, jpl DO_2D( ihls, ihls, ihls, ihls ) - !---------------------------------------- - ! zap ice energy and send it to the ocean - !---------------------------------------- - DO jk = 1, nlay_i - IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN - zhfx_res(ji,jj) = zhfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 - pe_i(ji,jj,jk,jl) = 0._wp - ENDIF - ENDDO - ! - DO jk = 1, nlay_s - IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN - zhfx_res(ji,jj) = zhfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 - pe_s(ji,jj,jk,jl) = 0._wp - ENDIF - ENDDO - ! - !----------------------------------------------------- - ! zap ice and snow volume, add water and salt to ocean - !----------------------------------------------------- IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN zwfx_res(ji,jj) = zwfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt pv_i (ji,jj,jl) = 0._wp @@ -672,7 +737,6 @@ CONTAINS pv_ip (ji,jj,jl) = 0._wp ENDIF END_2D - ! END DO ! ! record residual fluxes @@ -707,7 +771,6 @@ CONTAINS REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content !!------------------------------------------------------------------- - ! WHERE( pa_i (1:npti,:) < 0._wp ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0 WHERE( pv_i (1:npti,:) < 0._wp ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0 @@ -741,14 +804,15 @@ CONTAINS !!------------------------------------------------------------------- ! bv_i (:,:,:) = 0._wp + DO jl = 1, jpl + DO_3D( 0, 0, 0, 0, 1, nlay_i ) + IF( t_i(ji,jj,jk,jl) < rt0 - epsi10 ) THEN + bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rTmlt * sz_i(ji,jj,jk,jl) * r1_nlay_i / ( t_i(ji,jj,jk,jl) - rt0 ) + ENDIF + END_3D + ENDDO + ! DO_2D( 0, 0, 0, 0 ) - DO jl = 1, jpl - DO jk = 1, nlay_i - IF( t_i(ji,jj,jk,jl) < rt0 - epsi10 ) THEN - bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rTmlt * sz_i(ji,jj,jk,jl) * r1_nlay_i / ( t_i(ji,jj,jk,jl) - rt0 ) - ENDIF - ENDDO - ENDDO IF( vt_i(ji,jj) > epsi20 ) THEN bvm_i(ji,jj) = SUM( bv_i(ji,jj,:) * v_i(ji,jj,:) ) / vt_i(ji,jj) ELSE @@ -771,16 +835,19 @@ CONTAINS REAL(wp) :: ztmelts ! local scalar !!------------------------------------------------------------------- ! - DO ji = 1, npti - DO jk = 1, nlay_i ! Sea ice energy of melting - ztmelts = - rTmlt * sz_i_1d(ji,jk) + DO jk = 1, nlay_i ! Sea ice energy of melting + DO ji = 1, npti + ztmelts = - rTmlt * sz_i_1d(ji,jk) t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue ! (sometimes zdf scheme produces abnormally high temperatures) e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) & & + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) & & - rcp * ztmelts ) END DO - DO jk = 1, nlay_s ! Snow energy of melting + END DO + ! + DO jk = 1, nlay_s ! Snow energy of melting + DO ji = 1, npti e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) END DO END DO diff --git a/src/ICE/icewri.F90 b/src/ICE/icewri.F90 index 3cefb4d6be21080ae904ea5e85ff880022e2816a..d72fd7d807fd1174fd3b6c5a8506889531769dfd 100644 --- a/src/ICE/icewri.F90 +++ b/src/ICE/icewri.F90 @@ -72,15 +72,29 @@ CONTAINS ! tresholds for outputs DO_2D( 0, 0, 0, 0 ) - zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice - zmsk05(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.05_wp ) ) ! 1 if 5% ice , 0 if less - zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less - zmsksn(ji,jj) = MAX( 0._wp , SIGN( 1._wp , vt_s(ji,jj) - epsi06 ) ) ! 1 if snow , 0 if no snow - DO jl = 1, jpl - zmsk00l(ji,jj,jl) = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) - zmsksnl(ji,jj,jl) = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi06 ) ) - ENDDO + IF( at_i(ji,jj) >= epsi06 ) THEN ; zmsk00(ji,jj) = 1._wp ! 1 if ice , 0 if no ice + ELSE ; zmsk00(ji,jj) = 0._wp + ENDIF + IF( at_i(ji,jj) >= 0.05_wp ) THEN ; zmsk05(ji,jj) = 1._wp ! 1 if 5% ice , 0 if less + ELSE ; zmsk05(ji,jj) = 0._wp + ENDIF + IF( at_i(ji,jj) >= 0.15_wp ) THEN ; zmsk15(ji,jj) = 1._wp ! 1 if 15% ice, 0 if less + ELSE ; zmsk15(ji,jj) = 0._wp + ENDIF + IF( vt_s(ji,jj) >= epsi06 ) THEN ; zmsksn(ji,jj) = 1._wp ! 1 if snow , 0 if no snow + ELSE ; zmsksn(ji,jj) = 0._wp + ENDIF END_2D + DO jl = 1, jpl + DO_2D( 0, 0, 0, 0 ) + IF( a_i(ji,jj,jl) >= epsi06 ) THEN ; zmsk00l(ji,jj,jl) = 1._wp ! 1 if ice , 0 if no ice + ELSE ; zmsk00l(ji,jj,jl) = 0._wp + ENDIF + IF( v_s(ji,jj,jl) >= epsi06 ) THEN ; zmsksnl(ji,jj,jl) = 1._wp ! 1 if snow , 0 if no snow + ELSE ; zmsksnl(ji,jj,jl) = 0._wp + ENDIF + END_2D + ENDDO !----------------- ! Standard outputs diff --git a/src/OCE/SBC/sbc_oce.F90 b/src/OCE/SBC/sbc_oce.F90 index 5df541eb5e5a27d6c98d0f79a43531d9bd796dfe..d4eea48910fd30a7eba454c5bf6fcba1c69b6d46 100644 --- a/src/OCE/SBC/sbc_oce.F90 +++ b/src/OCE/SBC/sbc_oce.F90 @@ -119,7 +119,7 @@ MODULE sbc_oce REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp , emp_b !: freshwater budget: volume flux [Kg/m2/s] REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx , sfx_b !: salt flux [PSS.kg/m2/s] REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_tot !: total E-P over ocean and ice [Kg/m2/s] - REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fmmflx !: freshwater budget: freezing/melting [Kg/m2/s] + REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfice !: ice-ocean freshwater budget (>0 to the ocean) [Kg/m2/s] REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnf , rnf_b !: river runoff [Kg/m2/s] REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwficb !: iceberg melting [Kg/m2/s] !! @@ -201,7 +201,7 @@ CONTAINS & qns_tot(A2D(0)) , qsr_tot(A2D(0)) , qsr_hc(A2D(0),jpk) , qsr_hc_b(A2D(0),jpk) , STAT=ierr(5) ) ! ALLOCATE( sbc_tsc(A2D(0),jpts) , sbc_tsc_b(A2D(0),jpts) , & - & sfx (A2D(0)) , sfx_b(A2D(0)) , emp_tot(A2D(0)), fmmflx(A2D(0)), fwficb(A2D(0)), & + & sfx (A2D(0)) , sfx_b(A2D(0)) , emp_tot(A2D(0)), fwfice(A2D(0)), fwficb(A2D(0)), & & wndm(A2D(0)) , taum (A2D(0)) , STAT=ierr(6) ) ! ALLOCATE( tprecip(A2D(0)) , sprecip(A2D(0)) , & diff --git a/src/OCE/SBC/sbcice_cice.F90 b/src/OCE/SBC/sbcice_cice.F90 index 22833698c6af80ffa4c91547a163d6be2dff0331..be182ec281b19ac9d627ce9f44f5e57efd558493 100644 --- a/src/OCE/SBC/sbcice_cice.F90 +++ b/src/OCE/SBC/sbcice_cice.F90 @@ -564,7 +564,7 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2) WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0) sfx(:,:)=ztmp2(:,:)*1000.0 emp(:,:)=emp(:,:)-ztmp1(:,:) - fmmflx(:,:) = ztmp1(:,:) !!Joakim edit + fwfice(:,:) = -ztmp1(:,:) !!Joakim edit CALL lbc_lnk( 'sbcice_cice', emp , 'T', 1.0_wp, sfx , 'T', 1.0_wp ) diff --git a/src/OCE/SBC/sbcmod.F90 b/src/OCE/SBC/sbcmod.F90 index 58b695eba89a0a820ce49c059f907d02b316b9ab..1c7a3f574b3e52650a4bb65a78b8fb0d9f66eb68 100644 --- a/src/OCE/SBC/sbcmod.F90 +++ b/src/OCE/SBC/sbcmod.F90 @@ -228,7 +228,7 @@ CONTAINS ENDIF ! sfx (:,:) = 0._wp !* salt flux due to freezing/melting - fmmflx(:,:) = 0._wp !* freezing minus melting flux + fwfice(:,:) = 0._wp !* ice-ocean freshwater flux cloud_fra(:,:) = pp_cldf !* cloud fraction over sea ice (used in si3) taum(:,:) = 0._wp !* wind stress module (needed in GLS in case of reduced restart) @@ -614,7 +614,7 @@ CONTAINS CALL iom_put( "empmr" , emp(A2D(0))-rnf(A2D(0)) ) ! upward water flux CALL iom_put( "empbmr" , emp_b(A2D(0))-rnf(A2D(0)) ) ! before upward water flux (for ssh in offline ) CALL iom_put( "saltflx", sfx ) ! downward salt flux - CALL iom_put( "fmmflx" , fmmflx ) ! Freezing-melting water flux + CALL iom_put( "fwfice" , fwfice ) ! ice-ocean freshwater flux CALL iom_put( "qt" , qns+qsr ) ! total heat flux CALL iom_put( "qns" , qns ) ! solar heat flux CALL iom_put( "qsr" , qsr ) ! solar heat flux diff --git a/src/OFF/dtadyn.F90 b/src/OFF/dtadyn.F90 index 19fc6a1f39ab8271bc9a32b20214f98d95137bf7..d43cee979dddcf4157b391d84de9c77df6f4bea5 100644 --- a/src/OFF/dtadyn.F90 +++ b/src/OFF/dtadyn.F90 @@ -80,7 +80,7 @@ MODULE dtadyn INTEGER , SAVE :: jf_wnd ! index of wind speed INTEGER , SAVE :: jf_ice ! index of sea ice cover INTEGER , SAVE :: jf_rnf ! index of river runoff - INTEGER , SAVE :: jf_fmf ! index of downward salt flux + INTEGER , SAVE :: jf_fwf ! index of downward freshwater flux INTEGER , SAVE :: jf_ubl ! index of u-bbl coef INTEGER , SAVE :: jf_vbl ! index of v-bbl coef INTEGER , SAVE :: jf_div ! index of e3t @@ -136,13 +136,13 @@ CONTAINS ! IF( l_ldfslp .AND. .NOT.ln_c1d ) CALL dta_dyn_slp( kt, Kbb, Kmm ) ! Computation of slopes ! - ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask (:,:,:) ! temperature - ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask (:,:,:) ! salinity - wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * smask0(:,:) ! wind speed - needed for gas exchange - fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * smask0(:,:) ! downward salt flux (v3.5+) - fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask (:,:,1) ! Sea-ice fraction - qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * smask0(:,:) ! solar radiation - emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask (:,:,1) ! E-P + ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask (:,:,:) ! temperature + ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask (:,:,:) ! salinity + wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * smask0(:,:) ! wind speed - needed for gas exchange + fwfice(:,:) = - sf_dyn(jf_fwf)%fnow(:,:,1) * smask0(:,:) ! ice-ocean freshwater flux (>0 to the ocean) + fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask (:,:,1) ! Sea-ice fraction + qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * smask0(:,:) ! solar radiation + emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask (:,:,1) ! E-P IF( ln_dynrnf ) THEN rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask (:,:,1) ! rnf IF( ln_dynrnf_depth .AND. .NOT. ln_linssh ) CALL dta_dyn_rnf( Kmm ) @@ -232,14 +232,14 @@ CONTAINS TYPE(FLD_N), DIMENSION(jpfld) :: slf_d ! array of namelist informations on the fields to read TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp ! informations about the fields to be read TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt ! " " - TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf ! " " + TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fwf ! " " TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf ! " " TYPE(FLD_N) :: sn_div ! informations about the fields to be read !! NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, & & sn_uwd, sn_vwd, sn_wwd, sn_emp, & & sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr , & - & sn_wnd, sn_ice, sn_fmf, & + & sn_wnd, sn_ice, sn_fwf, & & sn_ubl, sn_vbl, sn_rnf, & & sn_empb, sn_div !!---------------------------------------------------------------------- @@ -263,13 +263,13 @@ CONTAINS ! jf_uwd = 1 ; jf_vwd = 2 ; jf_wwd = 3 ; jf_emp = 4 ; jf_avt = 5 jf_tem = 6 ; jf_sal = 7 ; jf_mld = 8 ; jf_qsr = 9 - jf_wnd = 10 ; jf_ice = 11 ; jf_fmf = 12 ; jfld = jf_fmf + jf_wnd = 10 ; jf_ice = 11 ; jf_fwf = 12 ; jfld = jf_fwf ! slf_d(jf_uwd) = sn_uwd ; slf_d(jf_vwd) = sn_vwd ; slf_d(jf_wwd) = sn_wwd slf_d(jf_emp) = sn_emp ; slf_d(jf_avt) = sn_avt slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal ; slf_d(jf_mld) = sn_mld slf_d(jf_qsr) = sn_qsr ; slf_d(jf_wnd) = sn_wnd ; slf_d(jf_ice) = sn_ice - slf_d(jf_fmf) = sn_fmf + slf_d(jf_fwf) = sn_fwf ! IF( .NOT.ln_linssh ) THEN jf_div = jfld + 1 ; jf_empb = jfld + 2 ; jfld = jf_empb @@ -313,7 +313,7 @@ CONTAINS IF( idimv == 3 ) THEN ; ipk = 1 ! 2D variable ELSE ; ipk = jpk ! 3D variable ENDIF - IF( ifpr == jf_mld .OR. ifpr == jf_qsr .OR. ifpr == jf_wnd .OR. ifpr == jf_fmf .OR. ifpr == jf_avt ) THEN + IF( ifpr == jf_mld .OR. ifpr == jf_qsr .OR. ifpr == jf_wnd .OR. ifpr == jf_fwf .OR. ifpr == jf_avt ) THEN ALLOCATE( sf_dyn(ifpr)%fnow(A2D(0),ipk) , STAT=ierr0 ) IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(A2D(0),ipk,2) , STAT=ierr1 ) ELSE diff --git a/src/TOP/PISCES/P4Z/p4zbc.F90 b/src/TOP/PISCES/P4Z/p4zbc.F90 index e4b20adbd59b360dfc4616924ab1d37023c6332c..c05ae9a10d4ad3891846c094efcf35d81dce094a 100644 --- a/src/TOP/PISCES/P4Z/p4zbc.F90 +++ b/src/TOP/PISCES/P4Z/p4zbc.F90 @@ -194,15 +194,15 @@ CONTAINS ! ------------------------------------------------------------------ DO_2D( 0, 0, 0, 0 ) zdep = rfact / e3t(ji,jj,1,Kmm) - zwflux = fmmflx(ji,jj) / 1000._wp - zironice = MAX( -0.99 * tr(ji,jj,1,jpfer,Kbb), -zwflux * icefeinput * zdep ) + zwflux = fwfice(ji,jj) / 1000._wp + zironice = MAX( -0.99 * tr(ji,jj,1,jpfer,Kbb), zwflux * icefeinput * zdep ) tr(ji,jj,1,jpfer,Krhs) = tr(ji,jj,1,jpfer,Krhs) + zironice END_2D ! IF( lk_iomput ) THEN ! iron flux from ice ALLOCATE( zw2d(A2D(0)) ) - zw2d(A2D(0)) = MAX( -0.99 * tr(A2D(0),1,jpfer,Kbb), (-1.*fmmflx(A2D(0))/1000._wp )*icefeinput*1.e+3*tmask(A2D(0),1)) + zw2d(A2D(0)) = MAX( -0.99 * tr(A2D(0),1,jpfer,Kbb), (fwfice(:,:)/1000._wp )*icefeinput*1.e+3*tmask(A2D(0),1)) CALL iom_put( "Ironice", zw2d ) DEALLOCATE( zw2d ) ENDIF diff --git a/src/TOP/TRP/trcsbc.F90 b/src/TOP/TRP/trcsbc.F90 index 88fd617f1efbeeb03ff099a56348c5dd3f6c1c30..42efe65bc125d0611295e7bfc0fdbec7ea05d174 100644 --- a/src/TOP/TRP/trcsbc.F90 +++ b/src/TOP/TRP/trcsbc.F90 @@ -51,12 +51,12 @@ CONTAINS !! * concentration/dilution effect: !! The surface freshwater flux modify the ocean volume !! and thus the concentration of a tracer as : - !! tr(Krhs) = tr(Krhs) + emp * tr(Kmm) / e3t_ + fmmflx * tri / e3t for k=1 + !! tr(Krhs) = tr(Krhs) + emp * tr(Kmm) / e3t_ + fwfice * tri / e3t for k=1 !! - tr(Kmm) , the concentration of tracer in the ocean !! - tri, the concentration of tracer in the sea-ice - !! - emp, the surface freshwater budget (evaporation minus precipitation + fmmflx) + !! - emp, the surface freshwater budget (evaporation minus precipitation + fwfice) !! given in kg/m2/s is divided by 1035 kg/m3 (density of ocean water) to obtain m/s. - !! - fmmflx, the flux asscociated to freezing-melting of sea-ice + !! - fwfice, the flux asscociated to freezing-melting of sea-ice !! In linear free surface case (ln_linssh=T), the volume of the !! ocean does not change with the water exchanges at the (air+ice)-sea !! @@ -132,7 +132,7 @@ CONTAINS ! DO jn = 1, jptra DO_2D( 0, 0, 0, 0 ) - sbc_trc(ji,jj,jn) = - fmmflx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) + sbc_trc(ji,jj,jn) = fwfice(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) END_2D END DO ! @@ -148,7 +148,7 @@ CONTAINS ! DO jn = 1, jptra DO_2D( 0, 0, 0, 0 ) - sbc_trc(ji,jj,jn) = - fmmflx(ji,jj) * r1_rho0 * trc_i(ji,jj,jn) + sbc_trc(ji,jj,jn) = fwfice(ji,jj) * r1_rho0 * trc_i(ji,jj,jn) END_2D END DO ! @@ -259,7 +259,7 @@ CONTAINS ! IF( .NOT.ln_linssh ) THEN !* only passive tracer fluxes associated with mass fluxes ! ! no passive tracer concentration modification due to ssh variation -!!st emp includes fmm see iceupdate.F90 +!!st emp includes fwfice see iceupdate.F90 !!not sure about trc_i case... (1) DO jn = 1, jptra DO_2D( 0, 0, 0, 0 ) !!st WHY 1 : exterior here ? @@ -292,12 +292,12 @@ CONTAINS END_2D END DO ! - CASE ( 0 ) ! Same concentration in sea ice and in the ocean fmm contribution to concentration/dilution effect has to be removed + CASE ( 0 ) ! Same concentration in sea ice and in the ocean fwfice contribution to concentration/dilution effect has to be removed ! DO jn = 1, jptra DO_2D( 0, 0, 0, 0 ) z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm) - ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + ( emp(ji,jj) - fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) + ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + ( emp(ji,jj) + fwfice(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) END_2D END DO ! @@ -307,13 +307,13 @@ CONTAINS DO_2D( 0, 0, 0, 0 ) z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm) ! tracer flux at the ice/ocean interface (tracer/m2/s) - zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice - ! ! only used in the levitating sea ice case + zftra = trc_i(ji,jj,jn) * fwfice(ji,jj) ! uptake of tracer in the sea ice + ! ! only used in the levitating sea ice case ! tracer flux only : add concentration dilution term in net tracer flux, no F-M in volume flux ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux ztfx = zftra ! net tracer flux ! - zdtra = r1_rho0 * ( ztfx + ( emp(ji,jj) - fmmflx(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) ) + zdtra = r1_rho0 * ( ztfx + ( emp(ji,jj) + fwfice(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) ) IF ( zdtra < 0. ) THEN zdtra = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc ) ! avoid negative concentrations to arise ENDIF @@ -333,7 +333,7 @@ CONTAINS DO jn = 1, jptra DO_2D( 0, 0, 0, 0 ) z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm) - ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) - fmmflx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) + ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + fwfice(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) END_2D END DO ! @@ -342,13 +342,13 @@ CONTAINS DO jn = 1, jptra DO_2D( 0, 0, 0, 0 ) ! tracer flux at the ice/ocean interface (tracer/m2/s) - zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice - ! ! only used in the levitating sea ice case + zftra = trc_i(ji,jj,jn) * fwfice(ji,jj) ! uptake of tracer in the sea ice + ! ! only used in the levitating sea ice case ! tracer flux only : add concentration dilution term in net tracer flux, no F-M in volume flux ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux ztfx = zftra ! net tracer flux ! - zdtra = r1_rho0 * ( ztfx - fmmflx(ji,jj) * ptr(ji,jj,1,jn,Kmm) ) + zdtra = r1_rho0 * ( ztfx + fwfice(ji,jj) * ptr(ji,jj,1,jn,Kmm) ) IF ( zdtra < 0. ) THEN zdtra = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc ) ! avoid negative concentrations to arise ENDIF diff --git a/src/TOP/oce_trc.F90 b/src/TOP/oce_trc.F90 index 7b3909107dd70d76dd1113d90389eede7fe1667c..5651a0e56a31651a4134988bae96bcf114a3d3a3 100644 --- a/src/TOP/oce_trc.F90 +++ b/src/TOP/oce_trc.F90 @@ -54,7 +54,7 @@ MODULE oce_trc USE sbc_oce , ONLY : qsr => qsr !: penetrative solar radiation (w m-2) USE sbc_oce , ONLY : emp => emp !: freshwater budget: volume flux [Kg/m2/s] USE sbc_oce , ONLY : emp_b => emp_b !: freshwater budget: volume flux [Kg/m2/s] - USE sbc_oce , ONLY : fmmflx => fmmflx !: freshwater budget: volume flux [Kg/m2/s] + USE sbc_oce , ONLY : fwfice => fwfice !: freshwater budget: volume flux [Kg/m2/s] USE sbc_oce , ONLY : rnf => rnf !: river runoff [Kg/m2/s] USE sbc_oce , ONLY : rnf_b => rnf_b !: river runoff at previus step [Kg/m2/s] USE sbc_oce , ONLY : ln_dm2dc => ln_dm2dc !: Diurnal Cycle