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 12949 – NEMO

Changeset 12949


Ignore:
Timestamp:
2020-05-18T23:58:40+02:00 (4 years ago)
Author:
clem
Message:

make sure ice and snow temperatures do not overshoot after advection (both with Prather or UMx), as well as ice salinity. Include the salinity bounds (rn_simin, rn_simax) inside the icethd_sal routine. There is no more problems with super cold temperatures now. Sette tests passed.

Location:
NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icecor.F90

    r11536 r12949  
    8181      DO jl = 1, jpl 
    8282         WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:) 
    83       END DO 
    84      
     83      END DO     
     84      !                             !----------------------------------------------------- 
     85      !                             !  Rebin categories with thickness out of bounds     ! 
     86      !                             !----------------------------------------------------- 
     87      IF ( jpl > 1 )   CALL ice_itd_reb( kt ) 
     88      ! 
    8589      !                             !----------------------------------------------------- 
    8690      IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        ! 
    … …  
    97101         END DO 
    98102      ENDIF 
    99       !                             !----------------------------------------------------- 
    100       !                             !  Rebin categories with thickness out of bounds     ! 
    101       !                             !----------------------------------------------------- 
    102       IF ( jpl > 1 )   CALL ice_itd_reb( kt ) 
    103  
    104103      !                             !----------------------------------------------------- 
    105104      CALL ice_var_zapsmall         !  Zap small values                                  ! 

  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_pra.F90

    r12832 r12949  
    8282      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    8383      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
    84       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il     ! melt pond lid thickness 
     84      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid thickness 
    8585      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8686      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
    … …  
    9292      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
    9393      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx 
    94       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max 
     94      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max, zs_i, zsi_max 
     95      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   ze_i, zei_max 
     96      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   ze_s, zes_max 
    9597      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea 
    9698      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
    … …  
    102104      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 
    103105      ! 
    104       ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
     106      ! --- Record max of the surrounding 9-pts (for call Hbig) --- ! 
     107      ! thickness and salinity 
     108      WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:) 
     109      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
     110      END WHERE 
    105111      DO jl = 1, jpl 
    106112         DO jj = 2, jpjm1 
    … …  
    118124                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    119125                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
     126               zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
     127                  &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
     128                  &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
     129                  &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
    120130            END DO 
    121131         END DO 
    122132      END DO 
    123       CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 
     133      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1., zsi_max, 'T', 1. ) 
     134      ! 
     135      ! enthalpies 
     136      DO jk = 1, nlay_i 
     137         WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:) 
     138         ELSEWHERE                      ; ze_i(:,:,jk,:) = 0._wp 
     139         END WHERE 
     140      END DO 
     141      DO jk = 1, nlay_s 
     142         WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:) 
     143         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
     144         END WHERE 
     145      END DO 
     146      DO jl = 1, jpl 
     147         DO jk = 1, nlay_i 
     148            DO jj = 2, jpjm1 
     149               DO ji = fs_2, fs_jpim1 
     150                  zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), & 
     151                     &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
     152                     &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
     153                     &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
     154               END DO 
     155            END DO 
     156         END DO 
     157      END DO 
     158      DO jl = 1, jpl 
     159         DO jk = 1, nlay_s 
     160            DO jj = 2, jpjm1 
     161               DO ji = fs_2, fs_jpim1 
     162                  zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), & 
     163                     &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
     164                     &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
     165                     &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
     166               END DO 
     167            END DO 
     168         END DO 
     169      END DO 
     170      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
     171      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
     172      ! 
    124173      ! 
    125174      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
    … …  
    283332         ! --- Make sure ice thickness is not too big --- ! 
    284333         !     (because ice thickness can be too large where ice concentration is very small) 
    285          CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     334         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 
     335            &            pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    286336         ! 
    287337         ! --- Ensure snow load is not too big --- ! 
    … …  
    635685 
    636686 
    637    SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     687   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 
     688      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    638689      !!------------------------------------------------------------------- 
    639690      !!                  ***  ROUTINE Hbig  *** 
    … …  
    649700      !! ** input   : Max thickness of the surrounding 9-points 
    650701      !!------------------------------------------------------------------- 
    651       REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
    652       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    653       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip 
     702      REAL(wp)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step 
     703      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts 
     704      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max 
     705      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max 
     706      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 
    654707      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
    655       ! 
    656       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    657       REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra 
     708      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     709      ! 
     710      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     711      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 
    658712      !!------------------------------------------------------------------- 
    659713      ! 
    … …  
    661715      ! 
    662716      DO jl = 1, jpl 
    663  
    664717         DO jj = 1, jpj 
    665718            DO ji = 1, jpi 
    … …  
    695748                  ENDIF            
    696749                  !                   
     750                  !                               ! -- check s_i -- ! 
     751                  ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 
     752                  zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 
     753                  IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     754                     zfra = psi_max(ji,jj,jl) / zsi 
     755                     sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 
     756                     psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 
     757                  ENDIF 
     758                  ! 
    697759               ENDIF 
    698760            END DO 
    699761         END DO 
    700762      END DO  
     763      ! 
     764      !                                           ! -- check e_i/v_i -- ! 
     765      DO jl = 1, jpl 
     766         DO jk = 1, nlay_i 
     767            DO jj = 1, jpj 
     768               DO ji = 1, jpi 
     769                  IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     770                     ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 
     771                     zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 
     772                     IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     773                        zfra = pei_max(ji,jj,jk,jl) / zei 
     774                        hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     775                        pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 
     776                     ENDIF 
     777                  ENDIF 
     778               END DO 
     779            END DO 
     780         END DO 
     781      END DO 
     782      !                                           ! -- check e_s/v_s -- ! 
     783      DO jl = 1, jpl 
     784         DO jk = 1, nlay_s 
     785            DO jj = 1, jpj 
     786               DO ji = 1, jpi 
     787                  IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 
     788                     ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 
     789                     zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 
     790                     IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     791                        zfra = pes_max(ji,jj,jk,jl) / zes 
     792                        hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     793                        pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 
     794                     ENDIF 
     795                  ENDIF 
     796               END DO 
     797            END DO 
     798         END DO 
     799      END DO 
    701800      ! 
    702801   END SUBROUTINE Hbig 

  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_umx.F90

    r12832 r12949  
    9393      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    9494      REAL(wp) ::   zdt, zvi_cen 
    95       REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication 
    96       REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
    97       REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2 
    98       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zu_cat, zv_cat 
    99       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups 
    100       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai , z1_aip, zhvar 
    101       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
     95      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
     96      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx, zcu_box, zcv_box 
     97      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
     98      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zu_cat, zv_cat 
     99      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zua_ho, zva_ho, zua_ups, zva_ups 
     100      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z1_ai , z1_aip, zhvar 
     101      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max, zs_i, zsi_max 
     102      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   ze_i, zei_max 
     103      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   ze_s, zes_max 
    102104      ! 
    103105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs  
    … …  
    106108      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    107109      ! 
    108       ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
     110      ! --- Record max of the surrounding 9-pts (for call Hbig) --- ! 
     111      ! thickness and salinity 
     112      WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:) 
     113      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
     114      END WHERE 
    109115      DO jl = 1, jpl 
    110116         DO jj = 2, jpjm1 
    … …  
    122128                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    123129                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    124             END DO 
    125          END DO 
    126       END DO 
    127       CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 
     130               zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
     131                  &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
     132                  &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
     133                  &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
     134            END DO 
     135         END DO 
     136      END DO 
     137      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1., zsi_max, 'T', 1. ) 
     138      ! 
     139      ! enthalpies 
     140      DO jk = 1, nlay_i 
     141         WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:) 
     142         ELSEWHERE                      ; ze_i(:,:,jk,:) = 0._wp 
     143         END WHERE 
     144      END DO 
     145      DO jk = 1, nlay_s 
     146         WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:) 
     147         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
     148         END WHERE 
     149      END DO 
     150      DO jl = 1, jpl 
     151         DO jk = 1, nlay_i 
     152            DO jj = 2, jpjm1 
     153               DO ji = fs_2, fs_jpim1 
     154                  zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), & 
     155                     &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
     156                     &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
     157                     &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
     158               END DO 
     159            END DO 
     160         END DO 
     161      END DO 
     162      DO jl = 1, jpl 
     163         DO jk = 1, nlay_s 
     164            DO jj = 2, jpjm1 
     165               DO ji = fs_2, fs_jpim1 
     166                  zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), & 
     167                     &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
     168                     &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
     169                     &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
     170               END DO 
     171            END DO 
     172         END DO 
     173      END DO 
     174      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
     175      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
    128176      ! 
    129177      ! 
    … …  
    362410         ! --- Make sure ice thickness is not too big --- ! 
    363411         !     (because ice thickness can be too large where ice concentration is very small) 
    364          CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     412         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 
     413            &            pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    365414         ! 
    366415         ! --- Ensure snow load is not too big --- ! 
    … …  
    15251574 
    15261575 
    1527    SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     1576   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 
     1577      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    15281578      !!------------------------------------------------------------------- 
    15291579      !!                  ***  ROUTINE Hbig  *** 
    … …  
    15391589      !! ** input   : Max thickness of the surrounding 9-points 
    15401590      !!------------------------------------------------------------------- 
    1541       REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
    1542       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    1543       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip 
     1591      REAL(wp)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step 
     1592      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts 
     1593      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max 
     1594      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max 
     1595      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 
    15441596      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
    1545       ! 
    1546       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    1547       REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra 
     1597      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     1598      ! 
     1599      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     1600      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 
    15481601      !!------------------------------------------------------------------- 
    15491602      ! 
    … …  
    15511604      ! 
    15521605      DO jl = 1, jpl 
    1553  
    15541606         DO jj = 1, jpj 
    15551607            DO ji = 1, jpi 
    … …  
    15851637                  ENDIF            
    15861638                  !                   
     1639                  !                               ! -- check s_i -- ! 
     1640                  ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 
     1641                  zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 
     1642                  IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1643                     zfra = psi_max(ji,jj,jl) / zsi 
     1644                     sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 
     1645                     psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 
     1646                  ENDIF 
     1647                  ! 
    15871648               ENDIF 
    15881649            END DO 
    15891650         END DO 
    15901651      END DO  
     1652      ! 
     1653      !                                           ! -- check e_i/v_i -- ! 
     1654      DO jl = 1, jpl 
     1655         DO jk = 1, nlay_i 
     1656            DO jj = 1, jpj 
     1657               DO ji = 1, jpi 
     1658                  IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     1659                     ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 
     1660                     zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 
     1661                     IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1662                        zfra = pei_max(ji,jj,jk,jl) / zei 
     1663                        hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1664                        pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 
     1665                     ENDIF 
     1666                  ENDIF 
     1667               END DO 
     1668            END DO 
     1669         END DO 
     1670      END DO 
     1671      !                                           ! -- check e_s/v_s -- ! 
     1672      DO jl = 1, jpl 
     1673         DO jk = 1, nlay_s 
     1674            DO jj = 1, jpj 
     1675               DO ji = 1, jpi 
     1676                  IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 
     1677                     ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 
     1678                     zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 
     1679                     IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1680                        zfra = pes_max(ji,jj,jk,jl) / zes 
     1681                        hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1682                        pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 
     1683                     ENDIF 
     1684                  ENDIF 
     1685               END DO 
     1686            END DO 
     1687         END DO 
     1688      END DO 
    15911689      ! 
    15921690   END SUBROUTINE Hbig 

  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_sal.F90

    r11536 r12949  
    5555      !!               -> nn_icesal = 3 -> Sice = S(z)   [multiyear ice] 
    5656      !!--------------------------------------------------------------------- 
    57       LOGICAL, INTENT(in) ::   ld_sal            ! gravity drainage and flushing or not  
     57      LOGICAL, INTENT(in) ::   ld_sal          ! gravity drainage and flushing or not  
    5858      ! 
    59       INTEGER  ::   ji, jk                       ! dummy loop indices  
    60       REAL(wp) ::   iflush, igravdr              ! local scalars 
    61       REAL(wp) ::   zs_sni, zs_i_gd, zs_i_fl, zs_i_si, zs_i_bg   ! local scalars 
     59      INTEGER  ::   ji                         ! dummy loop indices  
     60      REAL(wp) ::   zs_sni, zds                ! local scalars 
    6261      REAL(wp) ::   z1_time_gd, z1_time_fl 
    6362      !!--------------------------------------------------------------------- 
    … …  
    6867      CASE( 2 )       !  time varying salinity with linear profile  ! 
    6968         !            !---------------------------------------------! 
    70          z1_time_gd = 1._wp / rn_time_gd * rdt_ice 
    71          z1_time_fl = 1._wp / rn_time_fl * rdt_ice 
     69         z1_time_gd = rdt_ice / rn_time_gd 
     70         z1_time_fl = rdt_ice / rn_time_fl 
    7271         ! 
    7372         DO ji = 1, npti 
    7473            ! 
    75             !--------------------------------------------------------- 
    76             !  Update ice salinity from snow-ice and bottom growth 
    77             !--------------------------------------------------------- 
    7874            IF( h_i_1d(ji) > 0._wp ) THEN 
    79                zs_sni  = sss_1d(ji) * ( rhoi - rhos ) * r1_rhoi                     ! Salinity of snow ice 
    80                zs_i_si = ( zs_sni      - s_i_1d(ji) ) * dh_snowice(ji) / h_i_1d(ji) ! snow-ice     
    81                zs_i_bg = ( s_i_new(ji) - s_i_1d(ji) ) * dh_i_bog  (ji) / h_i_1d(ji) ! bottom growth 
    82                ! Update salinity (nb: salt flux already included in icethd_dh) 
    83                s_i_1d(ji) = s_i_1d(ji) + zs_i_bg + zs_i_si 
     75               ! 
     76               ! --- Update ice salinity from snow-ice and bottom growth --- ! 
     77               zs_sni = sss_1d(ji) * ( rhoi - rhos ) * r1_rhoi                           ! salinity of snow ice 
     78               zds    =       ( zs_sni      - s_i_1d(ji) ) * dh_snowice(ji) / h_i_1d(ji) ! snow-ice     
     79               zds    = zds + ( s_i_new(ji) - s_i_1d(ji) ) * dh_i_bog  (ji) / h_i_1d(ji) ! bottom growth 
     80               ! update salinity (nb: salt flux already included in icethd_dh) 
     81               s_i_1d(ji) = s_i_1d(ji) + zds 
     82               ! 
     83               ! --- Update ice salinity from brine drainage and flushing --- ! 
     84               IF( ld_sal ) THEN 
     85                  IF( t_su_1d(ji) >= rt0 ) THEN             ! flushing (summer time) 
     86                     zds = - MAX( s_i_1d(ji) - rn_sal_fl , 0._wp ) * z1_time_fl 
     87                  ELSEIF( t_su_1d(ji) <= t_bo_1d(ji) ) THEN ! gravity drainage 
     88                     zds = - MAX( s_i_1d(ji) - rn_sal_gd , 0._wp ) * z1_time_gd 
     89                  ELSE 
     90                     zds = 0._wp 
     91                  ENDIF 
     92                  ! update salinity 
     93                  s_i_1d(ji) = s_i_1d(ji) + zds 
     94                  ! salt flux 
     95                  sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoi * a_i_1d(ji) * h_i_1d(ji) * zds * r1_rdtice 
     96               ENDIF 
     97               ! 
     98               ! --- salinity must stay inbounds --- ! 
     99               zds =       MAX( 0._wp, rn_simin - s_i_1d(ji) ) ! > 0 if s_i < simin 
     100               zds = zds + MIN( 0._wp, rn_simax - s_i_1d(ji) ) ! < 0 if s_i > simax 
     101               ! update salinity 
     102               s_i_1d(ji) = s_i_1d(ji) + zds 
     103               ! salt flux 
     104               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * h_i_1d(ji) * zds * r1_rdtice 
     105               ! 
    84106            ENDIF 
    85107            ! 
    86             IF( ld_sal ) THEN 
    87                !--------------------------------------------------------- 
    88                !  Update ice salinity from brine drainage and flushing 
    89                !--------------------------------------------------------- 
    90                iflush   = MAX( 0._wp , SIGN( 1._wp , t_su_1d(ji) - rt0         ) )  ! =1 if summer  
    91                igravdr  = MAX( 0._wp , SIGN( 1._wp , t_bo_1d(ji) - t_su_1d(ji) ) )  ! =1 if t_su < t_bo 
    92  
    93                zs_i_gd = - igravdr * MAX( s_i_1d(ji) - rn_sal_gd , 0._wp ) * z1_time_gd  ! gravity drainage  
    94                zs_i_fl = - iflush  * MAX( s_i_1d(ji) - rn_sal_fl , 0._wp ) * z1_time_fl  ! flushing 
    95                 
    96                ! Update salinity    
    97                s_i_1d(ji) = s_i_1d(ji) + zs_i_fl + zs_i_gd 
    98                 
    99                ! Salt flux 
    100                sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoi * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_rdtice 
    101             ENDIF 
    102108         END DO 
    103109         ! 
Note: See TracChangeset for help on using the changeset viewer.