Ignore:
Timestamp:
2014-06-25T01:39:59+02:00 (6 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/limthd_sal.F90

    r4624 r4688  
    5353      ! 
    5454      INTEGER  ::   ji, jk     ! dummy loop indices  
    55       REAL(wp) ::   zsold, iflush, iaccrbo, igravdr, isnowic, i_ice_switch,  ztmelts   ! local scalars 
    56       REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars 
    57       REAL(wp), POINTER, DIMENSION(:) ::   ze_init, zhiold, zsiold 
     55      REAL(wp) ::   iflush, igravdr   ! local scalars 
    5856      !!--------------------------------------------------------------------- 
    5957 
    60       CALL wrk_alloc( jpij, ze_init, zhiold, zsiold ) 
    61  
     58      !--------------------------------------------------------- 
     59      !  0) Update ice salinity from snow-ice and bottom growth 
     60      !--------------------------------------------------------- 
     61      DO ji = kideb, kiut 
     62         sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji) 
     63      END DO 
     64  
    6265      !------------------------------------------------------------------------------| 
    6366      ! 1) Constant salinity, constant in time                                       | 
     
    7477      !  Module 2 : Constant salinity varying in time                                | 
    7578      !------------------------------------------------------------------------------| 
    76  
    7779      IF(  num_sal == 2  ) THEN 
    78  
    79          !--------------------------------- 
    80          ! Thickness at previous time step 
    81          !--------------------------------- 
    82          DO ji = kideb, kiut 
    83             zhiold(ji) = ht_i_b(ji) - dh_i_bott(ji) - dh_snowice(ji) - dh_i_surf(ji) 
    84             zsiold(ji) = sm_i_b(ji) 
    85          END DO 
    86  
    87          !--------------------- 
    88          ! Global heat content 
    89          !--------------------- 
    90          ze_init(:)  =  0._wp 
    91          DO jk = 1, nlay_i 
    92             DO ji = kideb, kiut 
    93                ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / REAL (nlay_i ) 
    94             END DO 
    95          END DO 
    9680 
    9781         DO ji = kideb, kiut 
     
    9983            ! Switches  
    10084            !---------- 
    101             iflush       =         MAX( 0._wp , SIGN( 1.0 , t_su_b(ji) - rtt )        )    ! =1 if summer  
    102             igravdr      =         MAX( 0._wp , SIGN( 1.0 , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo 
    103             iaccrbo      =         MAX( 0._wp , SIGN( 1.0 , dh_i_bott(ji) )           )    ! =1 if bottom accretion 
    104             i_ice_switch = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - ht_i_b(ji) + 1.e-2 ) ) 
    105             isnowic      = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - dh_snowice(ji) ) ) * i_ice_switch   ! =1 if snow ice formation 
     85            iflush  = MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt )        )    ! =1 if summer  
     86            igravdr = MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo 
    10687 
    10788            !--------------------- 
    10889            ! Salinity tendencies 
    10990            !--------------------- 
    110             !                                   ! drainage by gravity drainage 
     91            ! drainage by gravity drainage 
    11192            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice  
    112             !                                   ! drainage by flushing   
     93            ! drainage by flushing   
    11394            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice 
    11495 
     
    120101            sm_i_b(ji) = sm_i_b(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji) 
    121102 
    122             ! if no ice, salinity = 0.1 
    123             i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) ) 
    124             sm_i_b(ji)   = i_ice_switch * sm_i_b(ji) + s_i_min * ( 1._wp - i_ice_switch ) 
    125  
    126             !---------------------------- 
    127             ! Heat flux - brine drainage 
    128             !---------------------------- 
    129             fhbri_1d(ji) = 0._wp 
    130  
    131103            !---------------------------- 
    132104            ! Salt flux - brine drainage 
    133105            !---------------------------- 
    134             sfx_bri_1d(ji) = sfx_bri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) * r1_rdtice 
     106            sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoic * a_i_b(ji) * ht_i_b(ji) * ( dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji) ) * r1_rdtice 
    135107 
    136108         END DO 
     
    138110         ! Salinity profile 
    139111         CALL lim_var_salprof1d( kideb, kiut ) 
    140  
    141  
    142          ! Only necessary for conservation check since salinity is modified 
    143          !-------------------- 
    144          ! Temperature update 
    145          !-------------------- 
    146          DO jk = 1, nlay_i 
    147             DO ji = kideb, kiut 
    148                ztmelts    =  -tmut*s_i_b(ji,jk) + rtt 
    149                !Conversion q(S,T) -> T (second order equation) 
    150                zaaa         =  cpic 
    151                zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus 
    152                zccc         =  lfus * ( ztmelts - rtt ) 
    153                zdiscrim     =  SQRT(  MAX( zbbb*zbbb - 4.0*zaaa*zccc, 0._wp )  ) 
    154                t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa ) 
    155             END DO 
    156          END DO 
    157112         ! 
    158113      ENDIF  
     
    161116      !  Module 3 : Profile of salinity, constant in time                            | 
    162117      !------------------------------------------------------------------------------| 
    163  
    164118      IF(  num_sal == 3  )   CALL lim_var_salprof1d( kideb, kiut ) 
    165119 
    166       ! 
    167       CALL wrk_dealloc( jpij, ze_init, zhiold, zsiold ) 
    168120      ! 
    169121   END SUBROUTINE lim_thd_sal 
Note: See TracChangeset for help on using the changeset viewer.