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 13284 for NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_umx.F90 – NEMO

Ignore:
Timestamp:
2020-07-09T17:12:23+02:00 (4 years ago)
Author:
smasson
Message:

4.0-HEAD: merge 4.0-HEAD_r12713_clem_dan_fixcpl into 4.0-HEAD

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_umx.F90

    r12197 r13284  
    6060 
    6161   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  & 
    62       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     62      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    6363      !!---------------------------------------------------------------------- 
    6464      !!                  ***  ROUTINE ice_dyn_adv_umx  *** 
     
    8585      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond concentration 
    8686      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     87      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume 
    8788      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8889      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    9293      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    9394      REAL(wp) ::   zdt, zvi_cen 
    94       REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication 
    95       REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
    96       REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2 
    97       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zu_cat, zv_cat 
    98       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups 
    99       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai , z1_aip, zhvar 
    100       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 
    101104      ! 
    102105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs  
     
    105108      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    106109      ! 
    107       ! --- 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 
    108115      DO jl = 1, jpl 
    109116         DO jj = 2, jpjm1 
     
    121128                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    122129                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    123             END DO 
    124          END DO 
    125       END DO 
    126       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. ) 
    127176      ! 
    128177      ! 
     
    324373         ! 
    325374         !== melt ponds ==! 
    326          IF ( ln_pnd_H12 ) THEN 
     375         IF ( ln_pnd_LEV ) THEN 
    327376            ! concentration 
    328377            zamsk = 1._wp 
     
    334383            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
    335384               &                                      zhvar, pv_ip, zua_ups, zva_ups ) 
     385            ! lid 
     386            IF ( ln_pnd_lids ) THEN 
     387               zamsk = 0._wp 
     388               zhvar(:,:,:) = pv_il(:,:,:) * z1_aip(:,:,:) 
     389               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     390                  &                                      zhvar, pv_il, zua_ups, zva_ups ) 
     391            ENDIF 
    336392         ENDIF 
    337393         ! 
     
    350406         ! Remove negative values (conservation is ensured) 
    351407         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    352          CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     408         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    353409         ! 
    354410         ! --- Make sure ice thickness is not too big --- ! 
    355411         !     (because ice thickness can be too large where ice concentration is very small) 
    356          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 ) 
    357414         ! 
    358415         ! --- Ensure snow load is not too big --- ! 
     
    15171574 
    15181575 
    1519    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 ) 
    15201578      !!------------------------------------------------------------------- 
    15211579      !!                  ***  ROUTINE Hbig  *** 
     
    15311589      !! ** input   : Max thickness of the surrounding 9-points 
    15321590      !!------------------------------------------------------------------- 
    1533       REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
    1534       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    1535       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 
    15361596      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
    1537       ! 
    1538       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    1539       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 
    15401601      !!------------------------------------------------------------------- 
    15411602      ! 
     
    15431604      ! 
    15441605      DO jl = 1, jpl 
    1545  
    15461606         DO jj = 1, jpj 
    15471607            DO ji = 1, jpi 
     
    15501610                  !                               ! -- check h_ip -- ! 
    15511611                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    1552                   IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     1612                  IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    15531613                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    15541614                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     
    15771637                  ENDIF            
    15781638                  !                   
     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                  ! 
    15791648               ENDIF 
    15801649            END DO 
    15811650         END DO 
    15821651      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 
    15831689      ! 
    15841690   END SUBROUTINE Hbig 
Note: See TracChangeset for help on using the changeset viewer.