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 13472 for NEMO/trunk/src/ICE/icedyn_adv_umx.F90 – NEMO

Ignore:
Timestamp:
2020-09-16T15:05:19+02:00 (4 years ago)
Author:
smasson
Message:

trunk: commit changes from r4.0-HEAD from 13284 to 13449, see #2523

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_adv_umx.F90

    r13295 r13472  
    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_2D( 0, 0, 0, 0 ) 
     
    120127               &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    121128               &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    122          END_2D 
    123       END DO 
    124       CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp ) 
     129            zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
     130               &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
     131               &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
     132               &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
     133         END_2D 
     134      END DO 
     135      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
     136      ! 
     137      ! enthalpies 
     138      DO jk = 1, nlay_i 
     139         WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:) 
     140         ELSEWHERE                      ; ze_i(:,:,jk,:) = 0._wp 
     141         END WHERE 
     142      END DO 
     143      DO jk = 1, nlay_s 
     144         WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:) 
     145         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
     146         END WHERE 
     147      END DO 
     148      DO jl = 1, jpl 
     149         DO_3D( 0, 0, 0, 0, 1, nlay_i ) 
     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_3D 
     155      END DO 
     156      DO jl = 1, jpl 
     157         DO_3D( 0, 0, 0, 0, 1, nlay_s ) 
     158            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), & 
     159               &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
     160               &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
     161               &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
     162         END_3D 
     163      END DO 
     164      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
     165      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
    125166      ! 
    126167      ! 
     
    318359         ! 
    319360         !== melt ponds ==! 
    320          IF ( ln_pnd_H12 ) THEN 
     361         IF ( ln_pnd_LEV ) THEN 
    321362            ! concentration 
    322363            zamsk = 1._wp 
     
    328369            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
    329370               &                                      zhvar, pv_ip, zua_ups, zva_ups ) 
     371            ! lid 
     372            IF ( ln_pnd_lids ) THEN 
     373               zamsk = 0._wp 
     374               zhvar(:,:,:) = pv_il(:,:,:) * z1_aip(:,:,:) 
     375               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     376                  &                                      zhvar, pv_il, zua_ups, zva_ups ) 
     377            ENDIF 
    330378         ENDIF 
    331379         ! 
     
    342390         ! Remove negative values (conservation is ensured) 
    343391         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    344          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 ) 
     392         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 ) 
    345393         ! 
    346394         ! --- Make sure ice thickness is not too big --- ! 
    347395         !     (because ice thickness can be too large where ice concentration is very small) 
    348          CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     396         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 
     397            &            pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    349398         ! 
    350399         ! --- Ensure snow load is not too big --- ! 
     
    14091458 
    14101459 
    1411    SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     1460   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 
     1461      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    14121462      !!------------------------------------------------------------------- 
    14131463      !!                  ***  ROUTINE Hbig  *** 
     
    14231473      !! ** input   : Max thickness of the surrounding 9-points 
    14241474      !!------------------------------------------------------------------- 
    1425       REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
    1426       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    1427       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip 
     1475      REAL(wp)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step 
     1476      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts 
     1477      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max 
     1478      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max 
     1479      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 
    14281480      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
    1429       ! 
    1430       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    1431       REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra 
     1481      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     1482      ! 
     1483      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     1484      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 
    14321485      !!------------------------------------------------------------------- 
    14331486      ! 
     
    14351488      ! 
    14361489      DO jl = 1, jpl 
    1437  
    14381490         DO_2D( 1, 1, 1, 1 ) 
    14391491            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     
    14411493               !                               ! -- check h_ip -- ! 
    14421494               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    1443                IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     1495               IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    14441496                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    14451497                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     
    14681520               ENDIF            
    14691521               !                   
     1522               !                               ! -- check s_i -- ! 
     1523               ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 
     1524               zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 
     1525               IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1526                  zfra = psi_max(ji,jj,jl) / zsi 
     1527                  sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 
     1528                  psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 
     1529               ENDIF 
     1530               ! 
    14701531            ENDIF 
    14711532         END_2D 
    14721533      END DO  
     1534      ! 
     1535      !                                           ! -- check e_i/v_i -- ! 
     1536      DO jl = 1, jpl 
     1537         DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
     1538            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     1539               ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 
     1540               zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 
     1541               IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1542                  zfra = pei_max(ji,jj,jk,jl) / zei 
     1543                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1544                  pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 
     1545               ENDIF 
     1546            ENDIF 
     1547         END_3D 
     1548      END DO 
     1549      !                                           ! -- check e_s/v_s -- ! 
     1550      DO jl = 1, jpl 
     1551         DO_3D( 1, 1, 1, 1, 1, nlay_s ) 
     1552            IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 
     1553               ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 
     1554               zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 
     1555               IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1556                  zfra = pes_max(ji,jj,jk,jl) / zes 
     1557                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1558                  pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 
     1559               ENDIF 
     1560            ENDIF 
     1561         END_3D 
     1562      END DO 
    14731563      ! 
    14741564   END SUBROUTINE Hbig 
Note: See TracChangeset for help on using the changeset viewer.