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 for NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_pra.F90 – NEMO

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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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 
Note: See TracChangeset for help on using the changeset viewer.