New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9888 for NEMO/trunk/src/OCE – NEMO

Changeset 9888 for NEMO/trunk/src/OCE


Ignore:
Timestamp:
2018-07-06T12:33:01+02:00 (6 years ago)
Author:
clem
Message:

improve and debug BDY with sea ice. With this commit there should not be anymore problems in regional configurations (hopefully).

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/BDY/bdyice.F90

    r9885 r9888  
    139139         CALL lbc_bdy_lnk( h_s(:,:,jl), 'T', 1., ib_bdy ) 
    140140      ENDDO 
    141       ! retrieve at_i 
    142       at_i(:,:) = 0._wp 
    143       DO jl = 1, jpl 
    144          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    145       END DO 
    146141 
    147142      DO jl = 1, jpl 
     
    162157            !                                                             !      do not make state variables dependent on velocity 
    163158            ! 
    164             rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ii,ij) - 0.01 ) ) ! 0 if no ice 
    165             ! 
    166             ! concentration and thickness 
    167             a_i(ji,jj,jl) = a_i(ii,ij,jl) * rswitch 
    168             h_i(ji,jj,jl) = h_i(ii,ij,jl) * rswitch 
    169             h_s(ji,jj,jl) = h_s(ii,ij,jl) * rswitch 
    170             ! 
    171             ! Ice and snow volumes 
    172             v_i(ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) 
    173             v_s(ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) 
    174             ! 
    175             SELECT CASE( jpbound ) 
    176             ! 
    177             CASE( 0 )   ! velocity is inward 
    178                ! 
    179                ! Ice salinity, age, temperature 
    180                s_i (ji,jj,jl)   = rswitch * rn_ice_sal(ib_bdy)  + ( 1.0 - rswitch ) * rn_simin 
    181                oa_i(ji,jj,jl)   = rswitch * rn_ice_age(ib_bdy) * a_i(ji,jj,jl) 
    182                t_su(ji,jj,jl)   = rswitch * rn_ice_tem(ib_bdy)  + ( 1.0 - rswitch ) * rn_ice_tem(ib_bdy) 
     159            IF( a_i(ii,ij,jl) > 0._wp ) THEN   ! there is ice at the boundary 
     160               ! 
     161               a_i(ji,jj,jl) = a_i(ii,ij,jl) ! concentration 
     162               h_i(ji,jj,jl) = h_i(ii,ij,jl) ! thickness ice 
     163               h_s(ji,jj,jl) = h_s(ii,ij,jl) ! thickness snw 
     164               ! 
     165               SELECT CASE( jpbound ) 
     166                  ! 
     167               CASE( 0 )   ! velocity is inward 
     168                  ! 
     169                  oa_i(ji,jj,  jl) = rn_ice_age(ib_bdy) * a_i(ji,jj,jl) ! age 
     170                  a_ip(ji,jj,  jl) = 0._wp                              ! pond concentration 
     171                  v_ip(ji,jj,  jl) = 0._wp                              ! pond volume 
     172                  t_su(ji,jj,  jl) = rn_ice_tem(ib_bdy)                 ! temperature surface 
     173                  t_s (ji,jj,:,jl) = rn_ice_tem(ib_bdy)                 ! temperature snw 
     174                  t_i (ji,jj,:,jl) = rn_ice_tem(ib_bdy)                 ! temperature ice 
     175                  s_i (ji,jj,  jl) = rn_ice_sal(ib_bdy)                 ! salinity 
     176                  sz_i(ji,jj,:,jl) = rn_ice_sal(ib_bdy)                 ! salinity profile 
     177                  ! 
     178               CASE( 1 )   ! velocity is outward 
     179                  ! 
     180                  oa_i(ji,jj,  jl) = oa_i(ii,ij,  jl) ! age 
     181                  a_ip(ji,jj,  jl) = a_ip(ii,ij,  jl) ! pond concentration 
     182                  v_ip(ji,jj,  jl) = v_ip(ii,ij,  jl) ! pond volume 
     183                  t_su(ji,jj,  jl) = t_su(ii,ij,  jl) ! temperature surface 
     184                  t_s (ji,jj,:,jl) = t_s (ii,ij,:,jl) ! temperature snw 
     185                  t_i (ji,jj,:,jl) = t_i (ii,ij,:,jl) ! temperature ice 
     186                  s_i (ji,jj,  jl) = s_i (ii,ij,  jl) ! salinity 
     187                  sz_i(ji,jj,:,jl) = sz_i(ii,ij,:,jl) ! salinity profile 
     188                  ! 
     189               END SELECT 
     190               ! 
     191               IF( nn_icesal == 1 ) THEN     ! if constant salinity 
     192                  s_i (ji,jj  ,jl) = rn_icesal 
     193                  sz_i(ji,jj,:,jl) = rn_icesal 
     194               ENDIF 
     195               ! 
     196               ! global fields 
     197               v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)                       ! volume ice 
     198               v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)                       ! volume snw 
     199               sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 
    183200               DO jk = 1, nlay_s 
    184                   t_s(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0 
     201                  e_s(ji,jj,jk,jl) = rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )   ! enthalpy in J/m3 
     202                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s           ! enthalpy in J/m2 
     203               END DO                
     204               DO jk = 1, nlay_i 
     205                  ztmelts          = - tmut  * sz_i(ji,jj,jk,jl)              ! Melting temperature in C 
     206                  t_i(ji,jj,jk,jl) = MIN( t_i(ji,jj,jk,jl), ztmelts + rt0 )   ! Force t_i to be lower than melting point => likely conservation issue 
     207                  ! 
     208                  e_i(ji,jj,jk,jl) = rhoic * ( cpic * ( ztmelts - ( t_i(ji,jj,jk,jl) - rt0 ) )           &   ! enthalpy in J/m3 
     209                     &                       + lfus * ( 1._wp - ztmelts / ( t_i(ji,jj,jk,jl) - rt0 ) )   & 
     210                     &                       - rcp  *   ztmelts )                   
     211                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i                            ! enthalpy in J/m2 
    185212               END DO 
    186                DO jk = 1, nlay_i 
    187                   t_i (ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0  
    188                   sz_i(ji,jj,jk,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * rn_simin 
    189                END DO 
    190                ! 
    191                ! Ice ponds 
    192                a_ip(ji,jj,jl) = 0._wp 
    193                v_ip(ji,jj,jl) = 0._wp 
    194                ! 
    195             CASE( 1 )   ! velocity is outward 
    196                ! 
    197                ! Ice salinity, age, temperature 
    198                s_i (ji,jj,jl)   = rswitch * s_i (ii,ij,jl)  + ( 1.0 - rswitch ) * rn_simin 
    199                oa_i(ji,jj,jl)   = rswitch * oa_i(ii,ij,jl) 
    200                t_su(ji,jj,jl)   = rswitch * t_su(ii,ij,jl)  + ( 1.0 - rswitch ) * rt0 
    201                DO jk = 1, nlay_s 
    202                   t_s(ji,jj,jk,jl) = rswitch * t_s(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0 
    203                END DO 
    204                DO jk = 1, nlay_i 
    205                   t_i (ji,jj,jk,jl) = rswitch * t_i (ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0 
    206                   sz_i(ji,jj,jk,jl) = rswitch * sz_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rn_simin 
    207                END DO 
    208                ! 
    209                ! Ice ponds 
    210                a_ip(ji,jj,jl) = rswitch * a_ip(ii,ij,jl) 
    211                v_ip(ji,jj,jl) = rswitch * v_ip(ii,ij,jl) 
    212                ! 
    213             END SELECT 
    214             ! 
    215             IF( nn_icesal == 1 ) THEN     ! constant salinity : overwrite rn_icesal 
    216                s_i (ji,jj  ,jl) = rn_icesal 
    217                sz_i(ji,jj,:,jl) = rn_icesal 
     213               ! 
     214            ELSE   ! no ice at the boundary 
     215               ! 
     216               a_i (ji,jj,  jl) = 0._wp 
     217               h_i (ji,jj,  jl) = 0._wp 
     218               h_s (ji,jj,  jl) = 0._wp 
     219               oa_i(ji,jj,  jl) = 0._wp 
     220               a_ip(ji,jj,  jl) = 0._wp 
     221               v_ip(ji,jj,  jl) = 0._wp 
     222               t_su(ji,jj,  jl) = rt0 
     223               t_s (ji,jj,:,jl) = rt0 
     224               t_i (ji,jj,:,jl) = rt0  
     225                
     226               IF( nn_icesal == 1 ) THEN     ! if constant salinity 
     227                  s_i (ji,jj  ,jl) = rn_icesal 
     228                  sz_i(ji,jj,:,jl) = rn_icesal 
     229               ELSE                          ! if variable salinity 
     230                  s_i (ji,jj,jl)   = rn_simin 
     231                  sz_i(ji,jj,:,jl) = rn_simin 
     232               ENDIF 
     233               ! 
     234               ! global fields 
     235               v_i (ji,jj,  jl) = 0._wp 
     236               v_s (ji,jj,  jl) = 0._wp 
     237               sv_i(ji,jj,  jl) = 0._wp 
     238               e_s (ji,jj,:,jl) = 0._wp 
     239               e_i (ji,jj,:,jl) = 0._wp 
     240 
    218241            ENDIF 
    219             ! 
    220             ! contents 
    221             sv_i(ji,jj,jl)  = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 
    222             DO jk = 1, nlay_s 
    223                ! Snow energy of melting 
    224                e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 
    225                ! Multiply by volume, so that heat content in J/m2 
    226                e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
    227             END DO 
    228             DO jk = 1, nlay_i 
    229                ztmelts          = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K                   
    230                ! heat content per unit volume 
    231                e_i(ji,jj,jk,jl) = rswitch * rhoic * & 
    232                   (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    233                   +   lfus    * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 
    234                   - rcp      * ( ztmelts - rt0 ) ) 
    235                ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
    236                e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * h_i(ji,jj,jl) * r1_nlay_i 
    237             END DO 
    238             ! 
     242                         
    239243         END DO 
    240244         ! 
    241          CALL lbc_bdy_lnk( a_i(:,:,jl), 'T', 1., ib_bdy ) 
    242          CALL lbc_bdy_lnk( h_i(:,:,jl), 'T', 1., ib_bdy ) 
    243          CALL lbc_bdy_lnk( h_s(:,:,jl), 'T', 1., ib_bdy ) 
    244          CALL lbc_bdy_lnk( v_i(:,:,jl), 'T', 1., ib_bdy ) 
    245          CALL lbc_bdy_lnk( v_s(:,:,jl), 'T', 1., ib_bdy ) 
    246          ! 
     245         CALL lbc_bdy_lnk( a_i (:,:,jl), 'T', 1., ib_bdy ) 
     246         CALL lbc_bdy_lnk( h_i (:,:,jl), 'T', 1., ib_bdy ) 
     247         CALL lbc_bdy_lnk( h_s (:,:,jl), 'T', 1., ib_bdy ) 
     248         CALL lbc_bdy_lnk( oa_i(:,:,jl), 'T', 1., ib_bdy ) 
     249         CALL lbc_bdy_lnk( a_ip(:,:,jl), 'T', 1., ib_bdy ) 
     250         CALL lbc_bdy_lnk( v_ip(:,:,jl), 'T', 1., ib_bdy ) 
     251         CALL lbc_bdy_lnk( s_i (:,:,jl), 'T', 1., ib_bdy ) 
     252         CALL lbc_bdy_lnk( t_su(:,:,jl), 'T', 1., ib_bdy ) 
     253         CALL lbc_bdy_lnk( v_i (:,:,jl), 'T', 1., ib_bdy ) 
     254         CALL lbc_bdy_lnk( v_s (:,:,jl), 'T', 1., ib_bdy ) 
    247255         CALL lbc_bdy_lnk( sv_i(:,:,jl), 'T', 1., ib_bdy ) 
    248          CALL lbc_bdy_lnk(  s_i(:,:,jl), 'T', 1., ib_bdy ) 
    249          CALL lbc_bdy_lnk( oa_i(:,:,jl), 'T', 1., ib_bdy ) 
    250          CALL lbc_bdy_lnk( t_su(:,:,jl), 'T', 1., ib_bdy ) 
    251          DO jk = 1, nlay_s 
     256          DO jk = 1, nlay_s 
    252257            CALL lbc_bdy_lnk(t_s(:,:,jk,jl), 'T', 1., ib_bdy ) 
    253258            CALL lbc_bdy_lnk(e_s(:,:,jk,jl), 'T', 1., ib_bdy ) 
     
    258263         END DO 
    259264         ! 
    260          CALL lbc_bdy_lnk( a_ip(:,:,jl), 'T', 1., ib_bdy ) 
    261          CALL lbc_bdy_lnk( v_ip(:,:,jl), 'T', 1., ib_bdy ) 
    262          ! 
    263       END DO !jl 
     265      END DO ! jl 
    264266      !       
    265267   END SUBROUTINE bdy_ice_frs 
     
    309311                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice 
    310312                     !   
    311                      ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
     313                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then do not change u_ice) 
    312314                     u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 
    313315                        &            u_ice(ji-1,jj) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + & 
    314                         &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
     316                        &            u_ice(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
    315317                  ELSE                             ! everywhere else 
    316                      !u_ice(ji,jj) = u_oce(ji,jj) 
    317318                     u_ice(ji,jj) = 0._wp 
    318319                  ENDIF 
    319                   ! mask ice velocities 
    320                   rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01_wp ) ) ! 0 if no ice 
    321                   u_ice(ji,jj) = rswitch * u_ice(ji,jj) 
    322320                  ! 
    323321               END DO 
     
    336334                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice 
    337335                     !   
    338                      ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
     336                     ! v_ice = v_ice of the adjacent grid point except if this grid point is ice-free (then do not change v_ice) 
    339337                     v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 
    340338                        &            v_ice(ji,jj-1) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + & 
    341                         &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
     339                        &            v_ice(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
    342340                  ELSE                             ! everywhere else 
    343                      !v_ice(ji,jj) = v_oce(ji,jj) 
    344341                     v_ice(ji,jj) = 0._wp 
    345342                  ENDIF 
    346                   ! mask ice velocities 
    347                   rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01 ) ) ! 0 if no ice 
    348                   v_ice(ji,jj) = rswitch * v_ice(ji,jj) 
    349343                  ! 
    350344               END DO 
Note: See TracChangeset for help on using the changeset viewer.