Ignore:
Timestamp:
2014-06-25T01:39:59+02:00 (7 years ago)
Author:
clem
Message:

new version of LIM3 with perfect conservation of heat, see ticket #1352

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r4161 r4688  
    77   !!            3.0  ! 2004-06  (M. Vancoppenolle)   Energy Conservation  
    88   !!            4.0  ! 2011-02  (G. Madec)  add mpp considerations 
     9   !!             -   ! 2014-05  (C. Rousset) add lim_cons_hsm 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_lim3 
     
    1415   !!    lim_cons     :   checks whether energy, mass and salt are conserved  
    1516   !!---------------------------------------------------------------------- 
     17   USE phycst         ! physical constants 
    1618   USE par_ice        ! LIM-3 parameter 
    1719   USE ice            ! LIM-3 variables 
     
    2830   PUBLIC   lim_column_sum_energy 
    2931   PUBLIC   lim_cons_check 
     32   PUBLIC   lim_cons_hsm 
    3033 
    3134   !!---------------------------------------------------------------------- 
     
    151154   END SUBROUTINE lim_cons_check 
    152155 
     156 
     157   SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 
     158      !!------------------------------------------------------------------- 
     159      !!               ***  ROUTINE lim_cons_hsm *** 
     160      !! 
     161      !! ** Purpose : Test the conservation of heat, salt and mass for each routine 
     162      !! 
     163      !! ** Method  : 
     164      !!--------------------------------------------------------------------- 
     165      INTEGER         , INTENT(in)    :: icount      ! determine wether this is the beggining of the routine (0) or the end (1) 
     166      CHARACTER(len=*), INTENT(in)    :: cd_routine  ! name of the routine 
     167      REAL(wp)        , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
     168      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft 
     169      REAL(wp)                        :: zvmin, zamin, zamax  
     170 
     171      IF( icount == 0 ) THEN 
     172 
     173         zvi_b  = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
     174         zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     175         zei_b  = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
     176         zfw_b  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
     177            &                   wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) ) * area(:,:) * tms(:,:) ) 
     178         zfs_b  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
     179            &                   sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
     180         zft_b  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
     181            &                 - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
     182 
     183      ELSEIF( icount == 1 ) THEN 
     184 
     185         zfs  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
     186            &                 sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zfs_b 
     187         zfw  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
     188            &                 wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) ) * area(:,:) * tms(:,:) ) - zfw_b 
     189         zft  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
     190            &               - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zft_b 
     191  
     192         zvi  = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw  
     193         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 
     194         zei  =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zei_b * r1_rdtice + zft 
     195 
     196         zvmin = glob_min(v_i) 
     197         zamax = glob_max(SUM(a_i,dim=3)) 
     198         zamin = glob_min(a_i) 
     199        
     200         IF(lwp) THEN 
     201            IF ( ABS( zvi    ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday) 
     202            IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday) 
     203            IF ( ABS( zei    ) >  1.    ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (',cd_routine,') = ',(zei) 
     204            IF ( zvmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin) 
     205            IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+1.e-10 ) THEN 
     206                                          WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
     207            ENDIF 
     208            IF ( zamin <  0.            ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
     209         ENDIF 
     210 
     211      ENDIF 
     212 
     213   END SUBROUTINE lim_cons_hsm 
     214 
    153215#else 
    154216   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.