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 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 – NEMO

Ignore:
Timestamp:
2015-04-13T15:08:59+02:00 (9 years ago)
Author:
davestorkey
Message:

Merge in changes from trunk up to 5021.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r4688 r5208  
    2929   USE lib_mpp        ! MPP library 
    3030   USE wrk_nemo       ! work arrays 
     31   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3132   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3233   USE limthd_ent 
     
    3637 
    3738   PUBLIC lim_thd_lac     ! called by lim_thd 
    38  
    39    REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    40    REAL(wp) ::   epsi20 = 1.e-20_wp   ! 
    4139 
    4240   !!---------------------------------------------------------------------- 
     
    7169      !!             - Computation of variation of ice volume and mass 
    7270      !!             - Computation of frldb after lateral accretion and  
    73       !!               update ht_s_b, ht_i_b and tbif_1d(:,:)       
     71      !!               update ht_s_1d, ht_i_1d and tbif_1d(:,:)       
    7472      !!------------------------------------------------------------------------ 
    75       INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices 
    76       INTEGER ::   layer, nbpac     ! local integers  
    77       INTEGER ::   ii, ij, iter   !   -       - 
    78       REAL(wp)  ::   ztmelts, zdv, zfrazb, zweight, zindb, zinda, zde  ! local scalars 
     73      INTEGER ::   ji,jj,jk,jl      ! dummy loop indices 
     74      INTEGER ::   nbpac            ! local integers  
     75      INTEGER ::   ii, ij, iter     !   -       - 
     76      REAL(wp)  ::   ztmelts, zdv, zfrazb, zweight, zde                         ! local scalars 
    7977      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      - 
    8078      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      - 
     
    8987      REAL(wp) ::   zv_newfra 
    9088   
    91       INTEGER , POINTER, DIMENSION(:) ::   jcat      ! indexes of categories where new ice grows 
     89      INTEGER , POINTER, DIMENSION(:) ::   jcat        ! indexes of categories where new ice grows 
    9290      REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not 
    9391 
     
    10199      REAL(wp), POINTER, DIMENSION(:) ::   zda_res     ! residual area in case of excessive heat budget 
    102100      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_1d    ! total ice fraction     
    103       REAL(wp), POINTER, DIMENSION(:) ::   zat_i_lev   ! total ice fraction for level ice only (type 1)    
    104       REAL(wp), POINTER, DIMENSION(:) ::   zv_frazb   ! accretion of frazil ice at the ice bottom 
     101      REAL(wp), POINTER, DIMENSION(:) ::   zv_frazb    ! accretion of frazil ice at the ice bottom 
    105102      REAL(wp), POINTER, DIMENSION(:) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector) 
    106103 
    107       REAL(wp), POINTER, DIMENSION(:,:) ::   zv_old      ! old volume of ice in category jl 
    108       REAL(wp), POINTER, DIMENSION(:,:) ::   za_old      ! old area of ice in category jl 
    109       REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d     ! 1-D version of a_i 
    110       REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d     ! 1-D version of v_i 
    111       REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d    ! 1-D version of oa_i 
    112       REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d   ! 1-D version of smv_i 
     104      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_b      ! old volume of ice in category jl 
     105      REAL(wp), POINTER, DIMENSION(:,:) ::   za_b      ! old area of ice in category jl 
     106      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d   ! 1-D version of a_i 
     107      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d   ! 1-D version of v_i 
     108      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d  ! 1-D version of oa_i 
     109      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d ! 1-D version of smv_i 
    113110 
    114111      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d   !: 1-D version of e_i 
     
    119116      CALL wrk_alloc( jpij, jcat )   ! integer 
    120117      CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 
    121       CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d ) 
    122       CALL wrk_alloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
    123       CALL wrk_alloc( jpij,jkmax,jpl, ze_i_1d ) 
     118      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 
     119      CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
     120      CALL wrk_alloc( jpij,nlay_i+1,jpl, ze_i_1d ) 
    124121      CALL wrk_alloc( jpi,jpj, zvrel ) 
    125122 
     
    132129               DO ji = 1, jpi 
    133130                  !Energy of melting q(S,T) [J.m-3] 
    134                   zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes 
    135                   e_i(ji,jj,jk,jl) = zindb * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i ) 
     131                  rswitch          = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes 
     132                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) & 
     133                      &   / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i, wp ) 
    136134                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 
    137135               END DO 
     
    171169         zgamafr = 0.03 
    172170 
    173          DO jj = 1, jpj 
    174             DO ji = 1, jpi 
    175  
     171         DO jj = 2, jpj 
     172            DO ji = 2, jpi 
    176173               IF ( qlead(ji,jj) < 0._wp ) THEN 
    177174                  !------------- 
     
    189186                  ! Frazil ice velocity 
    190187                  !--------------------- 
    191                   zindb = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 
    192                   zvfrx = zindb * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 
    193                   zvfry = zindb * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 
     188                  rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 
     189                  zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 
     190                  zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 
    194191 
    195192                  !------------------- 
     
    197194                  !------------------- 
    198195                  ! C-grid ice velocity 
    199                   zindb = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  ) 
    200                   zvgx  = zindb * (  u_ice(ji-1,jj  ) * tmu(ji-1,jj  )    & 
    201                      &             + u_ice(ji,jj    ) * tmu(ji  ,jj  )  ) * 0.5_wp 
    202                   zvgy  = zindb * (  v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)    & 
    203                      &             + v_ice(ji,jj    ) * tmv(ji  ,jj  )  ) * 0.5_wp 
     196                  rswitch = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  ) 
     197                  zvgx    = rswitch * ( u_ice(ji-1,jj  ) * tmu(ji-1,jj  )  + u_ice(ji,jj) * tmu(ji,jj) ) * 0.5_wp 
     198                  zvgy    = rswitch * ( v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)  + v_ice(ji,jj) * tmv(ji,jj) ) * 0.5_wp 
    204199 
    205200                  !----------------------------------- 
     
    243238            END DO ! loop on ji ends 
    244239         END DO ! loop on jj ends 
     240      !  
     241      CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 
     242      CALL lbc_lnk( hicol(:,:), 'T', 1. ) 
    245243 
    246244      ENDIF ! End of computation of frazil ice collection thickness 
     
    255253      ! This occurs if open water energy budget is negative 
    256254      nbpac = 0 
     255      npac(:) = 0 
     256      ! 
    257257      DO jj = 1, jpj 
    258258         DO ji = 1, jpi 
     
    298298 
    299299         CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , jpi, jpj, npac(1:nbpac) ) 
    300          CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) ) 
     300         CALL tab_2d_1d( nbpac, t_bo_1d   (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) ) 
    301301         CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac)     , sfx_opw, jpi, jpj, npac(1:nbpac) ) 
    302302         CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) ) 
    303          CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) ) 
    304          CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) ) 
     303         CALL tab_2d_1d( nbpac, hicol_1d  (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) ) 
    305304         CALL tab_2d_1d( nbpac, zvrel_1d  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) ) 
    306305 
     
    315314         ! Keep old ice areas and volume in memory 
    316315         !----------------------------------------- 
    317          zv_old(:,:) = zv_i_1d(:,:)  
    318          za_old(:,:) = za_i_1d(:,:) 
    319  
     316         zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:)  
     317         za_b(1:nbpac,:) = za_i_1d(1:nbpac,:) 
    320318         !---------------------- 
    321319         ! Thickness of new ice 
     
    324322            zh_newice(ji) = hiccrit 
    325323         END DO 
    326          IF( fraz_swi == 1 ) zh_newice(:) = hicol_b(:) 
     324         IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 
    327325 
    328326         !---------------------- 
     
    331329         SELECT CASE ( num_sal ) 
    332330         CASE ( 1 )                    ! Sice = constant  
    333             zs_newice(:) = bulk_sal 
     331            zs_newice(1:nbpac) = bulk_sal 
    334332         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)] 
    335333            DO ji = 1, nbpac 
     
    339337            END DO 
    340338         CASE ( 3 )                    ! Sice = F(z) [multiyear ice] 
    341             zs_newice(:) =   2.3 
     339            zs_newice(1:nbpac) =   2.3 
    342340         END SELECT 
    343341 
     
    348346         DO ji = 1, nbpac 
    349347            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K) 
    350             ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             & 
    351                &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_b(ji) - rtt, -epsi10 ) )   & 
     348            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                             & 
     349               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) )   & 
    352350               &                       - rcp  *         ( ztmelts - rtt )  ) 
    353351         END DO ! ji 
     
    367365            zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg] 
    368366 
    369             zEw           = rcp * ( t_bo_b(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_b [J/kg] 
     367            zEw           = rcp * ( t_bo_1d(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_1d [J/kg] 
    370368                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)  
    371369                                                                    
     
    388386 
    389387            ! A fraction zfrazb of frazil ice is accreted at the ice bottom 
    390             zinda         = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 
    391             zfrazb        = zinda * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 
     388            rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 
     389            zfrazb        = rswitch * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 
    392390            zv_frazb(ji)  =         zfrazb   * zv_newice(ji) 
    393391            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 
     
    438436         DO ji = 1, nbpac 
    439437            jl = jcat(ji)                                                    ! categroy in which new ice is put 
    440             zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) ) )   ! 0 if old ice 
     438            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice 
    441439         END DO 
    442440 
     
    444442            DO ji = 1, nbpac 
    445443               jl = jcat(ji) 
    446                zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 
     444               rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 
    447445               ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                      & 
    448                   &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_old(ji,jl) )  & 
    449                   &        * zinda / MAX( zv_i_1d(ji,jl), epsi20 ) 
     446                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) )  & 
     447                  &        * rswitch / MAX( zv_i_1d(ji,jl), epsi20 ) 
    450448            END DO 
    451449         END DO 
     
    468466            ! new volumes including lateral/bottom accretion + residual 
    469467            DO ji = 1, nbpac 
    470                zinda          = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 
    471                zv_newfra      = zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 
    472                za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl)                
     468               rswitch        = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 
     469               zv_newfra      = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 
     470               za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl)                
    473471               zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 
    474  
    475472               ! for remapping 
    476473               h_i_old (ji,nlay_i+1) = zv_newfra 
    477474               qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 
    478475            ENDDO 
    479  
    480476            ! --- Ice enthalpy remapping --- ! 
    481             IF( zv_newfra > 0._wp ) THEN 
    482                CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) )  
    483             ENDIF 
    484  
     477            CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) )  
    485478         ENDDO 
    486479 
     
    490483         DO jl = 1, jpl 
    491484            DO ji = 1, nbpac 
    492                zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes 
    493                zoa_i_1d(ji,jl)  = za_old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb    
     485               rswitch          = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes 
     486               zoa_i_1d(ji,jl)  = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch    
    494487            END DO  
    495488         END DO    
     
    500493         DO jl = 1, jpl 
    501494            DO ji = 1, nbpac 
    502                zdv   = zv_i_1d(ji,jl) - zv_old(ji,jl) 
     495               zdv   = zv_i_1d(ji,jl) - zv_b(ji,jl) 
    503496               zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 
    504497            END DO 
     
    519512         CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 
    520513         CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 
    521          CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 
    522514 
    523515         CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) 
     
    534526               DO ji = 1, jpi 
    535527                  ! heat content in Joules 
    536                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) * unit_fac )  
     528                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac )  
    537529               END DO 
    538530            END DO 
     
    543535      CALL wrk_dealloc( jpij, jcat )   ! integer 
    544536      CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 
    545       CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d ) 
    546       CALL wrk_dealloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
    547       CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_1d ) 
     537      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 
     538      CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
     539      CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d ) 
    548540      CALL wrk_dealloc( jpi,jpj, zvrel ) 
    549541      ! 
Note: See TracChangeset for help on using the changeset viewer.