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 4161 for branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 – NEMO

Ignore:
Timestamp:
2013-11-07T11:01:27+01:00 (10 years ago)
Author:
cetlod
Message:

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r3625 r4161  
    4646 
    4747   !!---------------------------------------------------------------------- 
    48    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     48   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    4949   !! $Id$ 
    5050   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    7878      !!               update ht_s_b, ht_i_b and tbif_1d(:,:)       
    7979      !!------------------------------------------------------------------------ 
    80       INTEGER  ::   ji,jj,jk,jl,jm   ! dummy loop indices 
    81       INTEGER  ::   layer, nbpac     ! local integers  
    82       INTEGER  ::   zji, zjj, iter   !   -       - 
    83       REAL(wp) ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde   ! local scalars 
     80      INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices 
     81      INTEGER ::   layer, nbpac     ! local integers  
     82      INTEGER ::   ii, ij, iter   !   -       - 
     83      REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde  ! local scalars 
    8484      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      - 
    8585      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      - 
    86       REAL(wp) ::   zcoef                                                       !   -      - 
    8786      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness 
    8887      CHARACTER (len = 15) :: fieldid 
     
    160159               DO ji = 1, jpi 
    161160                  !Energy of melting q(S,T) [J.m-3] 
    162                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * nlay_i 
     161                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * REAL( nlay_i ) 
    163162                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) )  )   !0 if no ice and 1 if yes 
    164163                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 
     
    342341         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)] 
    343342            DO ji = 1, nbpac 
    344                zji =   MOD( npac(ji) - 1 , jpi ) + 1 
    345                zjj =      ( npac(ji) - 1 ) / jpi + 1 
    346                zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(zji,zjj)  ) 
     343               ii =   MOD( npac(ji) - 1 , jpi ) + 1 
     344               ij =      ( npac(ji) - 1 ) / jpi + 1 
     345               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij)  ) 
    347346            END DO 
    348347         CASE ( 3 )                    ! Sice = F(z) [multiyear ice] 
     
    389388         END DO 
    390389 
    391          !--------------------------------- 
    392          ! Salt flux due to new ice growth 
    393          !--------------------------------- 
    394          ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine) 
    395          DO ji = 1, nbpac 
    396             sfx_thd_1d(ji) = sfx_thd_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice 
    397             rdm_ice_1d(ji) = rdm_ice_1d(ji) +                 rhoic * zv_newice(ji) 
    398          END DO ! ji 
    399  
    400390         !------------------------------------ 
    401391         ! Diags for energy conservation test 
    402392         !------------------------------------ 
    403393         DO ji = 1, nbpac 
    404             zji = MOD( npac(ji) - 1 , jpi ) + 1 
    405             zjj =    ( npac(ji) - 1 ) / jpi + 1 
     394            ii = MOD( npac(ji) - 1 , jpi ) + 1 
     395            ij =    ( npac(ji) - 1 ) / jpi + 1 
    406396            ! 
    407             zde = ze_newice(ji) / unit_fac * area(zji,zjj) * zv_newice(ji) 
     397            zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji) 
    408398            ! 
    409             vt_i_init(zji,zjj) = vt_i_init(zji,zjj) + zv_newice(ji)             ! volume 
    410             et_i_init(zji,zjj) = et_i_init(zji,zjj) + zde                       ! Energy 
     399            vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji)             ! volume 
     400            et_i_init(ii,ij) = et_i_init(ii,ij) + zde                       ! Energy 
    411401 
    412402         END DO 
     
    419409         !----------------- 
    420410         DO ji = 1, nbpac 
    421             zji = MOD( npac(ji) - 1 , jpi ) + 1 
    422             zjj =    ( npac(ji) - 1 ) / jpi + 1 
     411            ii = MOD( npac(ji) - 1 , jpi ) + 1 
     412            ij =    ( npac(ji) - 1 ) / jpi + 1 
    423413            za_newice(ji) = zv_newice(ji) / zh_newice(ji) 
    424             diag_lat_gr(zji,zjj) = zv_newice(ji) * r1_rdtice 
     414            diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem 
    425415         END DO !ji 
    426416 
     
    441431         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
    442432         DO ji = 1, nbpac 
    443             IF ( za_newice(ji)  >  ( 1._wp - zat_i_ac(ji) ) ) THEN 
    444                zda_res(ji)   = za_newice(ji) - (1.0 - zat_i_ac(ji) ) 
     433            IF ( za_newice(ji) >  ( amax - zat_i_ac(ji) ) ) THEN 
     434               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_ac(ji) ) 
    445435               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji)  
    446436               za_newice(ji) = za_newice(ji) - zda_res  (ji) 
     
    473463         DO ji = 1, nbpac 
    474464            jl = zcatac(ji)                                                           ! categroy in which new ice is put 
    475             zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) ) )             ! zindb=1 if ice =0 otherwise 
     465            zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) )             ! zindb=1 if ice =0 otherwise 
    476466            zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb  ! old ice thickness 
    477467            zdhex    (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) )           ! difference in thickness 
    478             zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi11 ) )   ! ice totally new in jl category 
     468            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) )   ! ice totally new in jl category 
    479469         END DO 
    480470 
     
    482472            DO ji = 1, nbpac 
    483473               jl = zcatac(ji) 
    484                zqold   = ze_i_ac(ji,jk,jl)      ! [ J.m-3 ] 
    485                zalphai = MIN( zhice_old(ji,jl) *   jk       / nlay_i , zh_newice(ji) )   & 
    486                   &    - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i , zh_newice(ji) ) 
     474               zqold   = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 
     475               zalphai = MIN( zhice_old(ji,jl) * REAL( jk )     / REAL( nlay_i ), zh_newice(ji) )   & 
     476                  &    - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) ) 
    487477               ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji)                                     & 
    488                   + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / nlay_i   & 
     478                  + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / REAL( nlay_i )  & 
    489479                  + za_newice(ji)  * ze_newice(ji) * zalphai                                       & 
    490                   + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( zv_i_ac(ji,jl) / nlay_i ) 
     480                  + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) ) 
    491481            END DO 
    492482         END DO 
     
    513503            DO ji = 1, nbpac 
    514504               zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) ) 
    515                zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 ) 
     505               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi06 ) )  ! clem 
     506               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 ) 
    516507            END DO 
    517508         END DO 
     
    524515         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    525516            DO ji = 1, nbpac 
    526                zindb =  1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) ) )       ! zindb=1 if ice =0 otherwise 
     517               zindb =  1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) )       ! zindb=1 if ice =0 otherwise 
    527518               zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 
    528519               zdhicbot (ji,jl) = zdv_res(ji)    / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    & 
     
    536527            DO jk = 1, nlay_i 
    537528               DO ji = 1, nbpac 
    538                   zthick0(ji,jk,jl) =  zhice_old(ji,jl) / nlay_i 
     529                  zthick0(ji,jk,jl) =  zhice_old(ji,jl) / REAL( nlay_i ) 
    539530                  zqm0   (ji,jk,jl) =  ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl) 
    540531               END DO 
     
    555546               DO layer = 1, nlay_i + 1 
    556547                  DO ji = 1, nbpac 
    557                      zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) )  
     548                     zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )  
    558549                     ! Redistributing energy on the new grid 
    559                      zweight = MAX (  MIN( zhice_old(ji,jl) * layer , zdummy(ji,jl) * jk )   & 
    560                         &    - MAX( zhice_old(ji,jl) * ( layer - 1 ) , zdummy(ji,jl) * ( jk - 1 ) ) , 0._wp )   & 
    561                         &    /( MAX(nlay_i * zthick0(ji,layer,jl),epsi10) ) * zindb 
     550                     zweight = MAX (  MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) )   & 
     551                        &    - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp )   & 
     552                        &    /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb 
    562553                     ze_i_ac(ji,jk,jl) =  ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl)   
    563554                  END DO ! ji 
     
    569560            DO jk = 1, nlay_i 
    570561               DO ji = 1, nbpac 
    571                   zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )  
     562                  zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  
    572563                  ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl)   & 
    573                      &              / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * nlay_i * zindb 
     564                     &              / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb 
    574565               END DO 
    575566            END DO 
     
    581572         DO jl = 1, jpl 
    582573            DO ji = 1, nbpac 
    583                zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) )  ! 0 if no ice and 1 if yes 
     574               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes 
    584575               zoa_i_ac(ji,jl)  = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    
    585576            END DO  
     
    589580         ! Update salinity 
    590581         !----------------- 
    591          IF(  num_sal == 2  ) THEN      ! Sice = F(z,t) 
     582         !clem IF(  num_sal == 2  ) THEN 
    592583            DO jl = 1, jpl 
    593584               DO ji = 1, nbpac 
    594                   zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )   ! 0 if no ice and 1 if yes 
     585                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes 
    595586                  zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl) 
    596                   zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb 
     587                  zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif 
    597588               END DO 
    598589            END DO    
    599          ENDIF 
     590         !clem ENDIF 
     591 
     592         !-------------------------------- 
     593         ! Update mass/salt fluxes (clem) 
     594         !-------------------------------- 
     595         DO jl = 1, jpl 
     596            DO ji = 1, nbpac 
     597               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes 
     598               zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl) 
     599               rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb 
     600               sfx_thd_1d(ji)   =   sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb 
     601           END DO 
     602         END DO 
    600603 
    601604         !------------------------------------------------------------------------------! 
     
    606609            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 
    607610            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 
    608             IF (  num_sal == 2  )   & 
     611            !clem IF (  num_sal == 2  )   & 
    609612               CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 
    610613            DO jk = 1, nlay_i 
     
    622625      DO jl = 1, jpl 
    623626         DO jk = 1, nlay_i          ! heat content in 10^9 Joules 
    624             e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i / unit_fac  
     627            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac  
    625628         END DO 
    626629      END DO 
Note: See TracChangeset for help on using the changeset viewer.