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 5581 for branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90 – NEMO

Ignore:
Timestamp:
2015-07-10T13:28:53+02:00 (9 years ago)
Author:
timgraham
Message:

Merged head of trunk into branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    • Property svn:keywords set to Id
    r4688 r5581  
    55   !!====================================================================== 
    66   !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code 
    7    !!            3.6  !  2014-06  (C. Rousset)       Complete rewriting/cleaning 
     7   !!            3.5  !  2014-06  (C. Rousset)       Complete rewriting/cleaning 
    88   !!---------------------------------------------------------------------- 
    99#if defined key_lim3 
     
    1313   !!    lim_update1   : computes update of sea-ice global variables from trend terms 
    1414   !!---------------------------------------------------------------------- 
    15    USE limrhg          ! ice rheology 
    16  
    17    USE dom_oce 
    18    USE oce             ! dynamics and tracers variables 
    19    USE in_out_manager 
    2015   USE sbc_oce         ! Surface boundary condition: ocean fields 
    2116   USE sbc_ice         ! Surface boundary condition: ice fields 
    2217   USE dom_ice 
     18   USE dom_oce 
    2319   USE phycst          ! physical constants 
    2420   USE ice 
    25    USE limdyn 
    26    USE limtrp 
    27    USE limthd 
    28    USE limsbc 
    29    USE limdiahsb 
    30    USE limwri 
    31    USE limrst 
    3221   USE thd_ice         ! LIM thermodynamic sea-ice variables 
    33    USE par_ice 
    3422   USE limitd_th 
    35    USE limitd_me 
    3623   USE limvar 
    37    USE prtctl           ! Print control 
    38    USE lbclnk           ! lateral boundary condition - MPP exchanges 
    39    USE wrk_nemo         ! work arrays 
    40    USE lib_fortran     ! glob_sum 
    41    USE in_out_manager   ! I/O manager 
    42    USE iom              ! I/O manager 
    43    USE lib_mpp          ! MPP library 
     24   USE prtctl          ! Print control 
     25   USE wrk_nemo        ! work arrays 
    4426   USE timing          ! Timing 
    45    USE limcons        ! conservation tests 
     27   USE limcons         ! conservation tests 
     28   USE lib_mpp         ! MPP library 
     29   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     30   USE in_out_manager  ! I/O manager 
    4631 
    4732   IMPLICIT NONE 
    4833   PRIVATE 
    4934 
    50    PUBLIC   lim_update1   ! routine called by ice_step 
    51  
    52       REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       - 
    53           
     35   PUBLIC   lim_update1 
     36 
    5437   !! * Substitutions 
    5538#  include "vectopt_loop_substitute.h90" 
    5639   !!---------------------------------------------------------------------- 
    5740   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    58    !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $ 
     41   !! $Id$ 
    5942   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6043   !!---------------------------------------------------------------------- 
    6144CONTAINS 
    6245 
    63    SUBROUTINE lim_update1 
     46   SUBROUTINE lim_update1( kt ) 
    6447      !!------------------------------------------------------------------- 
    6548      !!               ***  ROUTINE lim_update1  *** 
     
    6952      !!                 
    7053      !!--------------------------------------------------------------------- 
    71       INTEGER  ::   ji, jj, jk, jl, jm    ! dummy loop indices 
    72       INTEGER  ::   jbnd1, jbnd2 
    73       INTEGER  ::   i_ice_switch 
     54      INTEGER, INTENT(in) ::   kt    ! number of iteration 
     55      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    7456      REAL(wp) ::   zsal 
    75       ! 
    76       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
     57      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    7758      !!------------------------------------------------------------------- 
    7859      IF( nn_timing == 1 )  CALL timing_start('limupdate1') 
     
    8061      IF( ln_limdyn ) THEN  
    8162 
     63      IF( kt == nit000 .AND. lwp ) THEN 
     64         WRITE(numout,*) ' lim_update1 '  
     65         WRITE(numout,*) ' ~~~~~~~~~~~ ' 
     66      ENDIF 
     67 
    8268      ! conservation test 
    8369      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    8470 
    85       !----------------- 
    86       ! zap small values 
    87       !----------------- 
    88       CALL lim_itd_me_zapsmall 
    89  
    90       CALL lim_var_glo2eqv 
    91       
    9271      !---------------------------------------------------- 
    93       ! Rebin categories with thickness out of bounds 
    94       !---------------------------------------------------- 
    95       DO jm = 1, jpm 
    96          jbnd1 = ice_cat_bounds(jm,1) 
    97          jbnd2 = ice_cat_bounds(jm,2) 
    98          IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 
    99       END DO 
    100  
     72      ! ice concentration should not exceed amax  
     73      !----------------------------------------------------- 
    10174      at_i(:,:) = 0._wp 
    10275      DO jl = 1, jpl 
     
    10477      END DO 
    10578 
    106       !---------------------------------------------------- 
    107       ! ice concentration should not exceed amax  
    108       !----------------------------------------------------- 
    10979      DO jl  = 1, jpl 
    11080         DO jj = 1, jpj 
    11181            DO ji = 1, jpi 
    112                IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    113                   a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) ) 
    114                   ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     82               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
     83                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     84                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
    11585               ENDIF 
    11686            END DO 
    11787         END DO 
    11888      END DO 
    119  
    120       at_i(:,:) = 0._wp 
    121       DO jl = 1, jpl 
    122          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    123       END DO 
    12489     
    125       ! -------------------------------------- 
    126       ! Final thickness distribution rebinning 
    127       ! -------------------------------------- 
    128       DO jm = 1, jpm 
    129          jbnd1 = ice_cat_bounds(jm,1) 
    130          jbnd2 = ice_cat_bounds(jm,2) 
    131          IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 
    132          IF (ice_ncat_types(jm) .EQ. 1 ) THEN 
    133          ENDIF 
    134       END DO 
    135  
    136       !----------------- 
    137       ! zap small values 
    138       !----------------- 
    139       CALL lim_itd_me_zapsmall 
    140  
    14190      !--------------------- 
    14291      ! Ice salinity bounds 
    14392      !--------------------- 
    144       IF (  num_sal == 2  ) THEN  
     93      IF (  nn_icesal == 2  ) THEN  
    14594         DO jl = 1, jpl 
    14695            DO jj = 1, jpj  
    14796               DO ji = 1, jpi 
    14897                  zsal            = smv_i(ji,jj,jl) 
    149                   smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl) 
    15098                  ! salinity stays in bounds 
    151                   i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 
    152                   smv_i(ji,jj,jl) = i_ice_switch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 
     99                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 
     100                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) 
    153101                  ! associated salt flux 
    154102                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
     
    158106      ENDIF 
    159107 
     108      !---------------------------------------------------- 
     109      ! Rebin categories with thickness out of bounds 
     110      !---------------------------------------------------- 
     111      IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 
     112 
     113      !----------------- 
     114      ! zap small values 
     115      !----------------- 
     116      CALL lim_var_zapsmall 
     117 
    160118      ! ------------------------------------------------- 
    161119      ! Diagnostics 
    162120      ! ------------------------------------------------- 
    163       d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:) 
    164       d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:) 
    165       d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:) 
    166       d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)   
    167       d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)    
    168       d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)   
    169       d_e_i_trp  (:,:,1:nlay_i,:) = e_i  (:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:) 
    170       d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
    171       d_smv_i_trp(:,:,:)   = 0._wp 
    172       IF(  num_sal == 2  ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     121      DO jl  = 1, jpl 
     122         afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 
     123      END DO 
     124 
     125      DO jj = 1, jpj 
     126         DO ji = 1, jpi             
     127            ! heat content variation (W.m-2) 
     128            diag_heat(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
     129               &                   SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
     130               &                 ) * r1_rdtice 
     131            ! salt, volume 
     132            diag_smvi(ji,jj) = SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 
     133            diag_vice(ji,jj) = SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice 
     134            diag_vsnw(ji,jj) = SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice 
     135         END DO 
     136      END DO 
    173137 
    174138      ! conservation test 
    175139      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    176140 
     141      ! ------------------------------------------------- 
     142      ! control prints 
     143      ! ------------------------------------------------- 
    177144      IF(ln_ctl) THEN   ! Control print 
    178145         CALL prt_ctl_info(' ') 
    179146         CALL prt_ctl_info(' - Cell values : ') 
    180147         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    181          CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update1  : cell area   :') 
     148         CALL prt_ctl(tab2d_1=e12t       , clinfo1=' lim_update1  : cell area   :') 
    182149         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update1  : at_i        :') 
    183150         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update1  : vt_i        :') 
     
    185152         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update1  : strength    :') 
    186153         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update1  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :') 
    187          CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update1  : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :') 
    188          CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update1  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :') 
     154         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update1  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :') 
    189155 
    190156         DO jl = 1, jpl 
     
    199165            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update1  : o_i         : ') 
    200166            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update1  : a_i         : ') 
    201             CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_a_i     : ') 
    202             CALL prt_ctl(tab2d_1=d_a_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_a_i_trp   : ') 
     167            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : a_i_b       : ') 
    203168            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update1  : v_i         : ') 
    204             CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_i     : ') 
    205             CALL prt_ctl(tab2d_1=d_v_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_i_trp   : ') 
     169            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_i_b       : ') 
    206170            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ') 
    207             CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_s     : ') 
    208             CALL prt_ctl(tab2d_1=d_v_s_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_s_trp   : ') 
    209             CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1        : ') 
    210             CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i1    : ') 
    211             CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : de_i1_trp   : ') 
    212             CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2        : ') 
    213             CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i2    : ') 
    214             CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : de_i2_trp   : ') 
     171            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_s_b       : ') 
     172            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_i1        : ') 
     173            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_i1_b      : ') 
     174            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)      , clinfo1= ' lim_update1  : e_i2        : ') 
     175            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)      , clinfo1= ' lim_update1  : e_i2_b      : ') 
    215176            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ') 
    216             CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update1  : old_e_snow  : ') 
    217             CALL prt_ctl(tab2d_1=d_e_s_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : d_e_s_trp   : ') 
     177            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow_b    : ') 
    218178            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update1  : smv_i       : ') 
    219             CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update1  : old_smv_i   : ') 
    220             CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl)        , clinfo1= ' lim_update1  : d_smv_i_trp : ') 
     179            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update1  : smv_i_b     : ') 
    221180            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update1  : oa_i        : ') 
    222             CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update1  : old_oa_i    : ') 
    223             CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ') 
     181            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update1  : oa_i_b      : ') 
    224182 
    225183            DO jk = 1, nlay_i 
Note: See TracChangeset for help on using the changeset viewer.