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 4649 for branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90 – NEMO

Ignore:
Timestamp:
2014-05-27T11:28:12+02:00 (10 years ago)
Author:
clem
Message:

finalizing LIM3 heat budget conservation + multiple minor bugs corrections

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r4634 r4649  
    3838   USE wrk_nemo         ! work arrays 
    3939   USE lib_fortran     ! glob_sum 
    40    ! Check budget (Rousset) 
    4140   USE in_out_manager   ! I/O manager 
    4241   USE iom              ! I/O manager 
    4342   USE lib_mpp          ! MPP library 
    4443   USE timing          ! Timing 
     44   USE limcons        ! conservation tests 
    4545 
    4646   IMPLICIT NONE 
     
    8282      REAL(wp) ::   zh, zdvres, zsal 
    8383      REAL(wp) ::   zat_i_old, ztmelts 
    84  
    85       REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset) 
    86       REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
     84      ! 
     85      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    8786      !!------------------------------------------------------------------- 
    8887      IF( nn_timing == 1 )  CALL timing_start('limupdate1') 
    8988 
    90       !------------------------------------------------------------------------------ 
    91       ! 1. Update of Global variables                                               | 
    92       !------------------------------------------------------------------------------ 
    93  
    94       !----------------- 
    95       !  Trend terms 
    96       !----------------- 
    97  
    98       ! ------------------------------- 
    99       !- check conservation (C Rousset) 
    100       IF (ln_limdiahsb) THEN 
    101          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
    102          zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    103          zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
    104          zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 
    105          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
    106          zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
    107       ENDIF 
    108       !- check conservation (C Rousset) 
    109       ! ------------------------------- 
    110  
     89      IF( ln_limdyn ) THEN  
     90 
     91      ! conservation test 
     92      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     93 
     94      !----------------- 
    11195      ! zap small values 
    11296      !----------------- 
     
    11599      CALL lim_var_glo2eqv 
    116100      
    117       at_i(:,:) = 0._wp 
    118       DO jl = 1, jpl 
    119          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    120       END DO 
    121  
    122101      !---------------------------------------------------- 
    123       ! 2.2) Rebin categories with thickness out of bounds 
     102      ! Rebin categories with thickness out of bounds 
    124103      !---------------------------------------------------- 
    125104      DO jm = 1, jpm 
     
    134113      END DO 
    135114 
    136       !--- 2.13 ice concentration should not exceed amax  
     115      !---------------------------------------------------- 
     116      ! ice concentration should not exceed amax  
    137117      !----------------------------------------------------- 
    138118      DO jl  = 1, jpl 
     
    147127      END DO 
    148128 
    149       at_i(:,:) = 0.0 
     129      at_i(:,:) = 0._wp 
    150130      DO jl = 1, jpl 
    151131         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    152132      END DO 
    153    
    154    
     133     
     134      ! -------------------------------------- 
    155135      ! Final thickness distribution rebinning 
    156136      ! -------------------------------------- 
     
    163143      END DO 
    164144 
    165  
     145      !----------------- 
    166146      ! zap small values 
    167147      !----------------- 
     
    169149 
    170150      !--------------------- 
    171       ! 2.11) Ice salinity 
     151      ! Ice salinity bounds 
    172152      !--------------------- 
    173153      IF (  num_sal == 2  ) THEN  
     
    179159                  ! salinity stays in bounds 
    180160                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 
    181                   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) ) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 
     161                  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) ) 
    182162                  ! associated salt flux 
    183163                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
    184                END DO ! ji 
    185             END DO ! jj 
    186          END DO !jl 
     164               END DO 
     165            END DO 
     166         END DO 
    187167      ENDIF 
    188  
    189       ! ------------------- 
    190       at_i(:,:) = a_i(:,:,1) 
    191       DO jl = 2, jpl 
    192          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    193       END DO 
    194    
    195168 
    196169      ! ------------------------------------------------- 
     
    206179      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
    207180      d_smv_i_trp(:,:,:)   = 0._wp 
    208       IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
    209  
    210       ! ------------------------------- 
    211       !- check conservation (C Rousset) 
    212       IF( ln_limdiahsb ) THEN 
    213          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    214          zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 
    215          zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 
    216   
    217          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw  
    218          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
    219          zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 
    220  
    221          zchk_vmin = glob_min(v_i) 
    222          zchk_amax = glob_max(SUM(a_i,dim=3)) 
    223          zchk_amin = glob_min(a_i) 
    224  
    225          IF(lwp) THEN 
    226             IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limupdate1) = ',(zchk_v_i * rday) 
    227             IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * rday) 
    228             IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limupdate1) = ',(zchk_e_i) 
    229             IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(zchk_vmin * 1.e-3) 
    230             IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate1) = ',zchk_amax 
    231             IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',zchk_amin 
    232          ENDIF 
    233        ENDIF 
    234       !- check conservation (C Rousset) 
    235       ! ------------------------------- 
     181      IF(  num_sal == 2  ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     182 
     183      ! conservation test 
     184      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    236185 
    237186      IF(ln_ctl) THEN   ! Control print 
     
    302251      ENDIF 
    303252    
     253      ENDIF ! ln_limdyn 
     254 
    304255      IF( nn_timing == 1 )  CALL timing_stop('limupdate1') 
    305256   END SUBROUTINE lim_update1 
Note: See TracChangeset for help on using the changeset viewer.