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 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 – NEMO

Ignore:
Timestamp:
2015-12-01T16:35:30+01:00 (8 years ago)
Author:
timgraham
Message:

Upgraded branch to r5518 of trunk (v3.6 stable revision)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r4333 r5965  
    1010   !!                 ! 2006-11 (X. Fettweis) Vectorized  
    1111   !!            3.0  ! 2008-03 (M. Vancoppenolle) Energy conservation and clean code 
    12    !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
     12   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
     13   !!             -   ! 2014-05 (C. Rousset) complete rewriting 
    1314   !!---------------------------------------------------------------------- 
    1415#if defined key_lim3 
     
    2223   USE domain         ! 
    2324   USE phycst         ! physical constants 
     25   USE sbc_oce        ! Surface boundary condition: ocean fields 
    2426   USE ice            ! LIM variables 
    25    USE par_ice        ! LIM parameters 
    2627   USE thd_ice        ! LIM thermodynamics 
    2728   USE limvar         ! LIM variables 
     
    3435   PRIVATE 
    3536 
    36    PUBLIC   lim_thd_ent         ! called by lim_thd 
    37  
    38    REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values 
    39    REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    40    REAL(wp) ::   zzero  = 0._wp      ! 
    41    REAL(wp) ::   zone   = 1._wp      ! 
     37   PUBLIC   lim_thd_ent         ! called by limthd and limthd_lac 
    4238 
    4339   !!---------------------------------------------------------------------- 
     
    4844CONTAINS 
    4945  
    50    SUBROUTINE lim_thd_ent( kideb, kiut, jl ) 
     46   SUBROUTINE lim_thd_ent( kideb, kiut, qnew ) 
    5147      !!------------------------------------------------------------------- 
    5248      !!               ***   ROUTINE lim_thd_ent  *** 
    5349      !! 
    5450      !! ** Purpose : 
    55       !!           This routine computes new vertical grids  
    56       !!           in the ice and in the snow, and consistently redistributes  
    57       !!           temperatures in the snow / ice.  
     51      !!           This routine computes new vertical grids in the ice,  
     52      !!           and consistently redistributes temperatures.  
    5853      !!           Redistribution is made so as to ensure to energy conservation 
    5954      !! 
     
    6156      !! ** Method  : linear conservative remapping 
    6257      !!            
    63       !! ** Steps : 1) Grid 
    64       !!            2) Switches 
    65       !!            3) Snow redistribution 
    66       !!            4) Ice enthalpy redistribution 
    67       !!            5) Ice salinity, recover temperature 
     58      !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses 
     59      !!            2) linear remapping on the new layers 
     60      !! 
     61      !! ------------ cum0(0)                        ------------- cum1(0) 
     62      !!                                    NEW      ------------- 
     63      !! ------------ cum0(1)               ==>      ------------- 
     64      !!     ...                                     ------------- 
     65      !! ------------                                ------------- 
     66      !! ------------ cum0(nlay_i+2)                 ------------- cum1(nlay_i) 
     67      !! 
    6868      !! 
    6969      !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 
    7070      !!------------------------------------------------------------------- 
    7171      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied 
    72       INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number 
    7372 
    74       INTEGER ::   ji,jk   !  dummy loop indices 
    75       INTEGER ::   ii, ij       ,   &  !  dummy indices 
    76          ntop0          ,   &  !  old layer top index 
    77          nbot1          ,   &  !  new layer bottom index 
    78          ntop1          ,   &  !  new layer top index 
    79          limsum         ,   &  !  temporary loop index 
    80          nlayi0,nlays0  ,   &  !  old number of layers 
    81          maxnbot0       ,   &  !  old layer bottom index 
    82          layer0, layer1        !  old/new layer indexes 
     73      REAL(wp), INTENT(inout), DIMENSION(:,:) :: qnew          ! new enthlapies (J.m-3, remapped) 
    8374 
    84  
    85       REAL(wp) :: & 
    86          ztmelts        ,   &  ! ice melting point 
    87          zqsnic         ,   &  ! enthalpy of snow ice layer 
    88          zhsnow         ,   &  ! temporary snow thickness variable 
    89          zswitch        ,   &  ! dummy switch argument 
    90          zfac1          ,   &  ! dummy factor 
    91          zfac2          ,   &  ! dummy factor 
    92          ztform         ,   &  !: bottom formation temperature 
    93          zaaa           ,   &  !: dummy factor 
    94          zbbb           ,   &  !: dummy factor 
    95          zccc           ,   &  !: dummy factor 
    96          zdiscrim              !: dummy factor 
    97  
    98       INTEGER, POINTER, DIMENSION(:) ::   snswi     !  snow switch 
    99       INTEGER, POINTER, DIMENSION(:) ::   nbot0     !  old layer bottom index 
    100       INTEGER, POINTER, DIMENSION(:) ::   icsuind   !  ice surface index 
    101       INTEGER, POINTER, DIMENSION(:) ::   icsuswi   !  ice surface switch 
    102       INTEGER, POINTER, DIMENSION(:) ::   icboind   !  ice bottom index 
    103       INTEGER, POINTER, DIMENSION(:) ::   icboswi   !  ice bottom switch 
    104       INTEGER, POINTER, DIMENSION(:) ::   snicind   !  snow ice index 
    105       INTEGER, POINTER, DIMENSION(:) ::   snicswi   !  snow ice switch 
    106       INTEGER, POINTER, DIMENSION(:) ::   snind     !  snow index 
     75      INTEGER  :: ji         !  dummy loop indices 
     76      INTEGER  :: jk0, jk1   !  old/new layer indices 
    10777      ! 
    108       REAL(wp), POINTER, DIMENSION(:) ::   zh_i   ! thickness of an ice layer 
    109       REAL(wp), POINTER, DIMENSION(:) ::   zh_s          ! thickness of a snow layer 
    110       REAL(wp), POINTER, DIMENSION(:) ::   zqsnow        ! enthalpy of the snow put in snow ice     
    111       REAL(wp), POINTER, DIMENSION(:) ::   zdeltah       ! temporary variable 
    112       REAL(wp), POINTER, DIMENSION(:) ::   zqti_in, zqts_in 
    113       REAL(wp), POINTER, DIMENSION(:) ::   zqti_fin, zqts_fin 
    114  
    115       REAL(wp), POINTER, DIMENSION(:,:) ::   zm0       !  old layer-system vertical cotes  
    116       REAL(wp), POINTER, DIMENSION(:,:) ::   qm0       !  old layer-system heat content  
    117       REAL(wp), POINTER, DIMENSION(:,:) ::   z_s       !  new snow system vertical cotes  
    118       REAL(wp), POINTER, DIMENSION(:,:) ::   z_i       !  new ice system vertical cotes  
    119       REAL(wp), POINTER, DIMENSION(:,:) ::   zthick0   !  old ice thickness  
    120       REAL(wp), POINTER, DIMENSION(:,:) ::   zhl0      ! old and new layer thicknesses  
    121       REAL(wp), POINTER, DIMENSION(:,:) ::   zrl01 
    122  
    123       REAL(wp) ::   zinda  
     78      REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces 
     79      REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces 
     80      REAL(wp), POINTER, DIMENSION(:)   :: zhnew               ! new layers thicknesses 
    12481      !!------------------------------------------------------------------- 
    12582 
    126       CALL wrk_alloc( jpij, snswi, nbot0, icsuind, icsuswi, icboind, icboswi, snicind, snicswi, snind )   ! integer 
    127       CALL wrk_alloc( jpij, zh_i, zh_s, zqsnow, zdeltah, zqti_in, zqts_in, zqti_fin, zqts_fin )           ! real 
    128       CALL wrk_alloc( jpij,jkmax+4, zm0, qm0, z_s, z_i, zthick0, zhl0, kjstart = 0 ) 
    129       CALL wrk_alloc( jkmax+4,jkmax+4, zrl01, kistart = 0, kjstart = 0 ) 
     83      CALL wrk_alloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 ) 
     84      CALL wrk_alloc( jpij, nlay_i+1, zqh_cum1, zh_cum1, kjstart = 0 ) 
     85      CALL wrk_alloc( jpij, zhnew ) 
    13086 
    131       zthick0(:,:) = 0._wp 
    132       zm0    (:,:) = 0._wp 
    133       qm0    (:,:) = 0._wp 
    134       zrl01  (:,:) = 0._wp 
    135       zhl0   (:,:) = 0._wp 
    136       z_i    (:,:) = 0._wp 
    137       z_s    (:,:) = 0._wp 
    138  
    139       ! 
    140       !------------------------------------------------------------------------------| 
    141       !  1) Grid                                                                     | 
    142       !------------------------------------------------------------------------------| 
    143       nlays0 = nlay_s 
    144       nlayi0 = nlay_i 
    145  
    146       DO ji = kideb, kiut 
    147          zh_i(ji) = old_ht_i_b(ji) / REAL( nlay_i )  
    148          zh_s(ji) = old_ht_s_b(ji) / REAL( nlay_s ) 
    149       END DO 
    150  
    151       ! 
    152       !------------------------------------------------------------------------------| 
    153       !  2) Switches                                                                 | 
    154       !------------------------------------------------------------------------------| 
    155       ! 2.1 snind(ji), snswi(ji) 
    156       ! snow surface behaviour : computation of snind(ji)-snswi(ji) 
    157       ! snind(ji) : index which equals  
    158       !   0 if snow is accumulating 
    159       !   1 if 1st layer is melting 
    160       !   2 if 2nd layer is melting ... 
    161       DO ji = kideb, kiut 
    162          snind  (ji) = 0 
    163          zdeltah(ji) = 0._wp 
    164       ENDDO !ji 
    165  
    166       DO jk = 1, nlays0 
     87      !-------------------------------------------------------------------------- 
     88      !  1) Cumulative integral of old enthalpy * thickness and layers interfaces 
     89      !-------------------------------------------------------------------------- 
     90      zqh_cum0(:,0:nlay_i+2) = 0._wp  
     91      zh_cum0 (:,0:nlay_i+2) = 0._wp 
     92      DO jk0 = 1, nlay_i+2 
    16793         DO ji = kideb, kiut 
    168             snind(ji)  = jk        *      NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)))) & 
    169                + snind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji))))) 
    170             zdeltah(ji)= zdeltah(ji) + zh_s(ji) 
    171          END DO ! ji 
    172       END DO ! jk 
    173  
    174       ! snswi(ji) : switch which value equals 1 if snow melts 
    175       !              0 if not 
    176       DO ji = kideb, kiut 
    177          snswi(ji)     = MAX(0,NINT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji))))) 
    178       END DO ! ji 
    179  
    180       ! 2.2 icsuind(ji), icsuswi(ji) 
    181       ! ice surface behaviour : computation of icsuind(ji)-icsuswi(ji) 
    182       ! icsuind(ji) : index which equals 
    183       !     0 if nothing happens at the surface 
    184       !     1 if first layer is melting 
    185       !     2 if 2nd layer is reached by melt ... 
    186       DO ji = kideb, kiut 
    187          icsuind(ji) = 0 
    188          zdeltah(ji) = 0._wp 
    189       END DO !ji 
    190       DO jk = 1, nlayi0 
    191          DO ji = kideb, kiut 
    192             icsuind(ji) = jk          *      NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)))) & 
    193                + icsuind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji))))) 
    194             zdeltah(ji) = zdeltah(ji) + zh_i(ji) 
    195          END DO ! ji 
    196       ENDDO !jk 
    197  
    198       ! icsuswi(ji) : switch which equals  
    199       !     1 if ice melts at the surface 
    200       !     0 if not 
    201       DO ji = kideb, kiut 
    202          icsuswi(ji)  = MAX(0,NINT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) ) 
     94            zqh_cum0(ji,jk0) = zqh_cum0(ji,jk0-1) + qh_i_old(ji,jk0-1) 
     95            zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + h_i_old (ji,jk0-1) 
     96         ENDDO 
    20397      ENDDO 
    20498 
    205       ! 2.3 icboind(ji), icboswi(ji) 
    206       ! ice bottom behaviour : computation of icboind(ji)-icboswi(ji) 
    207       ! icboind(ji) : index which equals 
    208       !     0 if accretion is on the way 
    209       !     1 if last layer has started to melt 
    210       !     2 if penultiem layer is melting ... and so on 
    211       !            N+1 if all layers melt and that snow transforms into ice 
    212       DO ji = kideb, kiut  
    213          icboind(ji) = 0 
    214          zdeltah(ji) = 0._wp 
    215       END DO 
    216       DO jk = nlayi0, 1, -1 
    217          DO ji = kideb, kiut 
    218             icboind(ji) = (nlayi0+1-jk) *      NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)))) & 
    219                &          + icboind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)))))  
    220             zdeltah(ji) = zdeltah(ji) + zh_i(ji) 
    221          END DO 
    222       END DO 
    223  
     99      !------------------------------------ 
     100      !  2) Interpolation on the new layers 
     101      !------------------------------------ 
     102      ! new layer thickesses 
    224103      DO ji = kideb, kiut 
    225          ! case of total ablation with remaining snow 
    226          IF ( ( ht_i_b(ji) .GT. epsi20 ) .AND. & 
    227             ( ht_i_b(ji) - dh_snowice(ji) .LT. epsi20 ) ) icboind(ji) = nlay_i + 1 
    228       END DO 
    229  
    230       ! icboswi(ji) : switch which equals  
    231       !     1 if ice accretion is on the way 
    232       !     0 if ablation is on the way 
    233       DO ji = kideb, kiut  
    234          icboswi(ji) = MAX(0,NINT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji))))) 
    235       END DO 
    236  
    237       ! 2.4 snicind(ji), snicswi(ji) 
    238       ! snow ice formation : calcul de snicind(ji)-snicswi(ji) 
    239       ! snicind(ji) : index which equals  
    240       !     0 if no snow-ice forms 
    241       !     1 if last layer of snow has started to melt 
    242       !     2 if penultiem layer ... 
    243       DO ji = kideb, kiut 
    244          snicind(ji) = 0 
    245          zdeltah(ji) = 0._wp 
    246       END DO 
    247       DO jk = nlays0, 1, -1 
    248          DO ji = kideb, kiut 
    249             snicind(ji) = (nlays0+1-jk) & 
    250                *      NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)))) + snicind(ji)   & 
    251                * (1 - NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji))))) 
    252             zdeltah(ji) = zdeltah(ji) + zh_s(ji) 
    253          END DO 
    254       END DO 
    255  
    256       ! snicswi(ji) : switch which equals  
    257       !     1 if snow-ice forms 
    258       !     0 if not 
    259       DO ji = kideb, kiut 
    260          snicswi(ji)   = MAX(0,NINT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji))))) 
     104         zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) * r1_nlay_i   
    261105      ENDDO 
    262106 
    263       ! 
    264       !------------------------------------------------------------------------------| 
    265       !  3) Snow redistribution                                                      | 
    266       !------------------------------------------------------------------------------| 
    267       ! 
    268       !------------- 
    269       ! Old profile 
    270       !------------- 
    271  
    272       ! by 'old', it is meant that layers coming from accretion are included,  
    273       ! and that interfacial layers which were partly melted are reduced  
    274  
    275       ! indexes of the vectors 
    276       !------------------------ 
    277       ntop0    =  1 
    278       maxnbot0 =  0 
    279  
    280       DO ji = kideb, kiut 
    281          nbot0(ji) =  nlays0  + 1 - snind(ji) + ( 1 - snicind(ji) ) * snicswi(ji) 
    282          ! cotes of the top of the layers 
    283          zm0(ji,0) =  0._wp 
    284          maxnbot0 =  MAX ( maxnbot0 , nbot0(ji) ) 
    285       END DO 
    286       IF( lk_mpp )   CALL mpp_max( maxnbot0, kcom=ncomm_ice ) 
    287  
    288       DO jk = 1, maxnbot0 
     107      ! new layers interfaces 
     108      zh_cum1(:,0:nlay_i) = 0._wp 
     109      DO jk1 = 1, nlay_i 
    289110         DO ji = kideb, kiut 
    290             !change 
    291             limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 
    292             limsum = MIN( limsum , nlay_s ) 
    293             zm0(ji,jk) =  dh_s_tot(ji) + zh_s(ji) * REAL( limsum ) 
    294          END DO 
    295       END DO 
    296  
    297       DO ji = kideb, kiut 
    298          zm0(ji,nbot0(ji)) =  dh_s_tot(ji) - REAL( snicswi(ji) ) * dh_snowice(ji) + zh_s(ji) * REAL( nlays0 ) 
    299          zm0(ji,1)         =  dh_s_tot(ji) * REAL( 1 - snswi(ji) ) + REAL( snswi(ji) ) * zm0(ji,1) 
    300       END DO 
    301  
    302       DO jk = ntop0, maxnbot0 
    303          DO ji = kideb, kiut 
    304             zthick0(ji,jk)  =  zm0(ji,jk) - zm0(ji,jk-1)            ! layer thickness 
    305          END DO 
    306       END DO 
    307  
    308       zqts_in(:) = 0._wp 
    309  
    310       DO ji = kideb, kiut         ! layer heat content 
    311          qm0    (ji,1) =  rhosn * (  cpic * ( rtt - REAL( 1 - snswi(ji) ) * tatm_ice_1d(ji)        & 
    312             &                                         - REAL( snswi(ji) ) * t_s_b      (ji,1)  )   & 
    313             &                      + lfus  ) * zthick0(ji,1) 
    314          zqts_in(ji)   =  zqts_in(ji) + qm0(ji,1)  
    315       END DO 
    316  
    317       DO jk = 2, maxnbot0 
    318          DO ji = kideb, kiut 
    319             limsum      = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 
    320             limsum      = MIN( limsum , nlay_s ) 
    321             qm0(ji,jk)  = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk) 
    322             zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, - ht_s_b(ji) ) ) 
    323             zqts_in(ji) = zqts_in(ji) + REAL( 1 - snswi(ji) ) * qm0(ji,jk) * zswitch 
    324          END DO ! jk 
    325       END DO ! ji 
    326  
    327       !------------------------------------------------ 
    328       ! Energy given by the snow in snow-ice formation 
    329       !------------------------------------------------ 
    330       ! zqsnow, enthalpy of the flooded snow 
    331       DO ji = kideb, kiut 
    332          zqsnow (ji) =  rhosn * lfus 
    333          zdeltah(ji) =  0._wp 
    334       END DO 
    335  
    336       DO jk =  nlays0, 1, -1 
    337          DO ji = kideb, kiut 
    338             zhsnow =  MAX( 0._wp , dh_snowice(ji)-zdeltah(ji) ) 
    339             zqsnow (ji) =  zqsnow (ji) + rhosn*cpic*(rtt-t_s_b(ji,jk)) 
    340             zdeltah(ji) =  zdeltah(ji) + zh_s(ji) 
    341          END DO 
    342       END DO 
    343  
    344       DO ji = kideb, kiut 
    345          zqsnow(ji) = zqsnow(ji) * dh_snowice(ji) 
    346       END DO 
    347  
    348       !------------------ 
    349       ! new snow profile 
    350       !------------------ 
    351  
    352       !-------------- 
    353       ! Vector index    
    354       !-------------- 
    355       ntop1 =  1 
    356       nbot1 =  nlay_s 
    357  
    358       !------------------- 
    359       ! Layer coordinates  
    360       !------------------- 
    361       DO ji = kideb, kiut 
    362          zh_s(ji)  = ht_s_b(ji) / REAL( nlay_s ) 
    363          z_s(ji,0) =  0._wp 
     111            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji) 
     112         ENDDO 
    364113      ENDDO 
    365114 
    366       DO jk = 1, nlay_s 
     115      zqh_cum1(:,0:nlay_i) = 0._wp  
     116      ! new cumulative q*h => linear interpolation 
     117      DO jk0 = 1, nlay_i+1 
     118         DO jk1 = 1, nlay_i-1 
     119            DO ji = kideb, kiut 
     120               IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN 
     121                  zqh_cum1(ji,jk1) = ( zqh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1  ) ) +  & 
     122                     &                 zqh_cum0(ji,jk0  ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) )  & 
     123                     &             / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) ) 
     124               ENDIF 
     125            ENDDO 
     126         ENDDO 
     127      ENDDO 
     128      ! to ensure that total heat content is strictly conserved, set: 
     129      zqh_cum1(:,nlay_i) = zqh_cum0(:,nlay_i+2)  
     130 
     131      ! new enthalpies 
     132      DO jk1 = 1, nlay_i 
    367133         DO ji = kideb, kiut 
    368             z_s(ji,jk) =  zh_s(ji) * REAL( jk ) 
    369          END DO 
    370       END DO 
    371  
    372       !----------------- 
    373       ! Layer thickness 
    374       !----------------- 
    375       DO layer0 = ntop0, maxnbot0 
    376          DO ji = kideb, kiut 
    377             zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) 
    378          END DO 
    379       END DO 
    380  
    381       DO layer1 = ntop1, nbot1 
    382          DO ji = kideb, kiut 
    383             q_s_b(ji,layer1) = 0._wp 
    384          END DO 
    385       END DO 
    386  
    387       !---------------- 
    388       ! Weight factors 
    389       !---------------- 
    390       DO layer0 = ntop0, maxnbot0 
    391          DO layer1 = ntop1, nbot1 
    392             DO ji = kideb, kiut 
    393                zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) 
    394                zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1))   & 
    395                   &                 - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10))  
    396                q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0)   & 
    397                   &                                * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 
    398             END DO 
    399          END DO 
    400       END DO 
    401  
    402       ! Heat conservation 
    403       zqts_fin(:) = 0._wp 
    404       DO jk = 1, nlay_s 
    405          DO ji = kideb, kiut 
    406             zqts_fin(ji) = zqts_fin(ji) + q_s_b(ji,jk) 
    407          END DO 
    408       END DO 
    409  
    410       IF ( con_i .AND. jiindex_1d > 0 ) THEN 
    411          DO ji = kideb, kiut 
    412             IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN 
    413                ii                 = MOD( npb(ji) - 1, jpi ) + 1 
    414                ij                 = ( npb(ji) - 1 ) / jpi + 1 
    415                WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 
    416                WRITE(numout,*) ' ji, jj   : ', ii, ij 
    417                WRITE(numout,*) ' ht_s_b   : ', ht_s_b(ji) 
    418                WRITE(numout,*) ' zqts_in  : ', zqts_in (ji) * r1_rdtice 
    419                WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) * r1_rdtice 
    420                WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji) 
    421                WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji) 
    422                WRITE(numout,*) ' snswi    : ', snswi(ji) 
    423             ENDIF 
    424          END DO 
    425       ENDIF 
    426  
    427       !--------------------- 
    428       ! Recover heat content 
    429       !--------------------- 
    430       DO jk = 1, nlay_s 
    431          DO ji = kideb, kiut 
    432             zinda = MAX( 0._wp, SIGN( 1._wp , zh_s(ji) - epsi10 ) )         
    433             q_s_b(ji,jk) = zinda * q_s_b(ji,jk) / MAX( zh_s(ji) , epsi10 ) 
    434          END DO !ji 
    435       END DO !jk   
    436  
    437       !--------------------- 
    438       ! Recover temperature 
    439       !--------------------- 
    440       zfac1 = 1. / ( rhosn * cpic ) 
    441       zfac2 = lfus / cpic   
    442       DO jk = 1, nlay_s 
    443          DO ji = kideb, kiut 
    444             zswitch = MAX ( 0.0 , SIGN ( 1.0, - ht_s_b(ji) ) ) 
    445             t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 ) 
    446          END DO 
    447       END DO 
    448       ! 
    449       !------------------------------------------------------------------------------| 
    450       !  4) Ice redistribution                                                       | 
    451       !------------------------------------------------------------------------------| 
    452       ! 
    453       !------------- 
    454       ! OLD PROFILE  
    455       !------------- 
    456  
    457       !---------------- 
    458       ! Vector indexes 
    459       !---------------- 
    460       ntop0    =  1 
    461       maxnbot0 =  0 
    462  
    463       DO ji = kideb, kiut 
    464          ! reference number of the bottommost layer 
    465          nbot0(ji) =  MAX( 1 ,  MIN( nlayi0 + ( 1 - icboind(ji) ) +        & 
    466             &                           ( 1 - icsuind(ji) ) * icsuswi(ji) + snicswi(ji) , nlay_i + 2 ) ) 
    467          ! maximum reference number of the bottommost layer over all domain 
    468          maxnbot0 =  MAX( maxnbot0 , nbot0(ji) ) 
    469       END DO 
    470  
    471       !------------------------- 
    472       ! Cotes of old ice layers 
    473       !------------------------- 
    474       zm0(:,0) =  0._wp 
    475  
    476       DO jk = 1, maxnbot0 
    477          DO ji = kideb, kiut 
    478             ! jk goes from 1 to nbot0 
    479             ! the ice layer number goes from 1 to nlay_i 
    480             ! limsum is the real ice layer number corresponding to present jk 
    481             limsum    =  ( (icsuswi(ji)*(icsuind(ji)+jk-1) + &  
    482                (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) 
    483             zm0(ji,jk)=  REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) & 
    484                +  REAL(limsum) * zh_i(ji) 
    485          END DO 
    486       END DO 
    487  
    488       DO ji = kideb, kiut 
    489          zm0(ji,nbot0(ji)) =  REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) + dh_i_bott(ji) & 
    490             +  zh_i(ji) * REAL(nlayi0) 
    491          zm0(ji,1)         =  REAL(snicswi(ji))*dh_snowice(ji) + REAL(1-snicswi(ji))*zm0(ji,1) 
    492       END DO 
    493  
    494       !----------------------------- 
    495       ! Thickness of old ice layers 
    496       !----------------------------- 
    497       DO jk = ntop0, maxnbot0 
    498          DO ji = kideb, kiut 
    499             zthick0(ji,jk) =  zm0(ji,jk) - zm0(ji,jk-1) 
    500          END DO 
    501       END DO 
    502  
    503       !--------------------------- 
    504       ! Inner layers heat content 
    505       !--------------------------- 
    506       qm0(:,:) =  0.0 
    507       zqti_in(:) = 0.0 
    508  
    509       DO jk = ntop0, maxnbot0 
    510          DO ji = kideb, kiut 
    511             limsum =  MAX(1,MIN(snicswi(ji)*(jk-1) + icsuswi(ji)*(jk-1+icsuind(ji)) + & 
    512                (1-icsuswi(ji))*(1-snicswi(ji))*jk,nlay_i)) 
    513             ztmelts = -tmut * s_i_b(ji,limsum) + rtt 
    514             qm0(ji,jk) = rhoic * ( cpic * (ztmelts-t_i_b(ji,limsum)) + lfus * ( 1.0-(ztmelts-rtt)/ & 
    515                MIN((t_i_b(ji,limsum)-rtt),-epsi20) ) - rcp*(ztmelts-rtt) ) & 
    516                * zthick0(ji,jk) 
    517          END DO 
    518       END DO 
    519  
    520       !---------------------------- 
    521       ! Bottom layers heat content 
    522       !---------------------------- 
    523       DO ji = kideb, kiut         
    524          ztmelts    = REAL( 1 - icboswi(ji) ) * (-tmut * s_i_b  (ji,nlayi0) )   &   ! case of melting ice 
    525             &       +     REAL( icboswi(ji) ) * (-tmut * s_i_new(ji)        )   &   ! case of forming ice 
    526             &       + rtt                                                         ! in Kelvin 
    527  
    528          ! bottom formation temperature 
    529          ztform = t_i_b(ji,nlay_i) 
    530          IF(  num_sal == 2  )   ztform = t_bo_b(ji) 
    531          qm0(ji,nbot0(ji)) = REAL( 1 - icboswi(ji) )*qm0(ji,nbot0(ji))             &   ! case of melting ice 
    532             &              + REAL( icboswi(ji) ) * rhoic * ( cpic*(ztmelts-ztform)       &   ! case of forming ice 
    533             + lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) )      &  
    534             - rcp*(ztmelts-rtt) ) * zthick0(ji,nbot0(ji)  ) 
    535       END DO 
    536  
    537       !----------------------------- 
    538       ! Snow ice layer heat content 
    539       !----------------------------- 
    540       DO ji = kideb, kiut 
    541          ! energy of the flooding seawater 
    542          zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * & 
    543             (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive 
    544          ! Heat conservation diagnostic 
    545          qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic  
    546  
    547          qldif_1d(ji)   = qldif_1d(ji) + zqsnic * a_i_b(ji) 
    548  
    549          ! enthalpy of the newly formed snow-ice layer 
    550          ! = enthalpy of snow + enthalpy of frozen water 
    551          zqsnic         =  zqsnow(ji) + zqsnic 
    552          qm0(ji,1)      =  REAL(snicswi(ji)) * zqsnic + REAL( 1 - snicswi(ji) ) * qm0(ji,1) 
    553  
    554       END DO ! ji 
    555  
    556       DO jk = ntop0, maxnbot0 
    557          DO ji = kideb, kiut 
    558             ! Heat conservation 
    559             zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi10) ) & 
    560                &                                   * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) 
    561          END DO 
    562       END DO 
    563  
    564       !------------- 
    565       ! NEW PROFILE 
    566       !------------- 
    567  
    568       !--------------- 
    569       ! Vectors index 
    570       !--------------- 
    571       ntop1 =  1  
    572       nbot1 =  nlay_i 
    573  
    574       !------------------ 
    575       ! Layers thickness  
    576       !------------------ 
    577       DO ji = kideb, kiut 
    578          zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 
     134            rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) )  
     135            qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 
     136         ENDDO 
    579137      ENDDO 
    580138 
    581       !------------- 
    582       ! Layer cotes       
    583       !------------- 
    584       z_i(:,0) =  0._wp 
    585       DO jk = 1, nlay_i 
    586          DO ji = kideb, kiut 
    587             z_i(ji,jk) =  zh_i(ji) * jk 
    588          END DO 
     139      ! --- diag error on heat remapping --- ! 
     140      ! comment: if input h_i_old and qh_i_old are already multiplied by a_i (as in limthd_lac),  
     141      ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0 
     142      DO ji = kideb, kiut 
     143         hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice *  & 
     144            &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( qh_i_old(ji,0:nlay_i+1) ) )  
    589145      END DO 
    590  
    591       !--thicknesses of the layers 
    592       DO layer0 = ntop0, maxnbot0 
    593          DO ji = kideb, kiut 
    594             zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1)   ! thicknesses of the layers 
    595          END DO 
    596       END DO 
    597  
    598       !------------------------ 
    599       ! Weights for relayering 
    600       !------------------------ 
    601       q_i_b(:,:) = 0._wp 
    602       DO layer0 = ntop0, maxnbot0 
    603          DO layer1 = ntop1, nbot1 
    604             DO ji = kideb, kiut 
    605                zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) 
    606                zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) & 
    607                   - MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10)) 
    608                q_i_b(ji,layer1) = q_i_b(ji,layer1) &  
    609                   + zrl01(layer1,layer0)*qm0(ji,layer0) & 
    610                   * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi10)) & 
    611                   * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 
    612             END DO 
    613          END DO 
    614       END DO 
    615  
    616       !------------------------- 
    617       ! Heat conservation check 
    618       !------------------------- 
    619       zqti_fin(:) = 0._wp 
    620       DO jk = 1, nlay_i 
    621          DO ji = kideb, kiut 
    622             zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk) 
    623          END DO 
    624       END DO 
     146       
    625147      ! 
    626       IF ( con_i .AND. jiindex_1d > 0 ) THEN 
    627          DO ji = kideb, kiut 
    628             IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN 
    629                ii                 = MOD( npb(ji) - 1, jpi ) + 1 
    630                ij                 = ( npb(ji) - 1 ) / jpi + 1 
    631                WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 
    632                WRITE(numout,*) ' ji, jj   : ', ii, ij 
    633                WRITE(numout,*) ' ht_i_b   : ', ht_i_b(ji) 
    634                WRITE(numout,*) ' zqti_in  : ', zqti_in (ji) * r1_rdtice 
    635                WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice 
    636                WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) 
    637                WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) 
    638                WRITE(numout,*) ' dh_snowice:', dh_snowice(ji) 
    639                WRITE(numout,*) ' icsuswi  : ', icsuswi(ji) 
    640                WRITE(numout,*) ' icboswi  : ', icboswi(ji) 
    641                WRITE(numout,*) ' snicswi  : ', snicswi(ji) 
    642             ENDIF 
    643          END DO 
    644       ENDIF 
    645  
    646       !---------------------- 
    647       ! Recover heat content  
    648       !---------------------- 
    649       DO jk = 1, nlay_i 
    650          DO ji = kideb, kiut 
    651             zinda = MAX( 0._wp, SIGN( 1._wp , zh_i(ji) - epsi10 ) ) 
    652             q_i_b(ji,jk) = zinda * q_i_b(ji,jk) / MAX( zh_i(ji) , epsi10 ) 
    653          END DO !ji 
    654       END DO !jk   
    655  
    656       ! Heat conservation 
    657       zqti_fin(:) = 0.0 
    658       DO jk = 1, nlay_i 
    659          DO ji = kideb, kiut 
    660             zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk) * zh_i(ji) 
    661          END DO 
    662       END DO 
    663  
    664       ! 
    665       !------------------------------------------------------------------------------| 
    666       !  5) Update salinity and recover temperature                                  | 
    667       !------------------------------------------------------------------------------| 
    668       ! 
    669       ! Update salinity (basal entrapment, snow ice formation) 
    670       DO ji = kideb, kiut 
    671          sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji) 
    672       END DO !ji 
    673  
    674       ! Recover temperature 
    675       DO jk = 1, nlay_i 
    676          DO ji = kideb, kiut 
    677             ztmelts    =  -tmut*s_i_b(ji,jk) + rtt 
    678             !Conversion q(S,T) -> T (second order equation) 
    679             zaaa         =  cpic 
    680             zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus 
    681             zccc         =  lfus * ( ztmelts - rtt ) 
    682             zdiscrim     =  SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) ) 
    683             t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa ) 
    684          END DO !ji 
    685  
    686       END DO !jk 
    687       ! 
    688       CALL wrk_dealloc( jpij, snswi, nbot0, icsuind, icsuswi, icboind, icboswi, snicind, snicswi, snind )   ! integer 
    689       CALL wrk_dealloc( jpij, zh_i, zh_s, zqsnow, zdeltah, zqti_in, zqts_in, zqti_fin, zqts_fin )           ! real 
    690       CALL wrk_dealloc( jpij,jkmax+4, zm0, qm0, z_s, z_i, zthick0, zhl0, kjstart = 0 ) 
    691       CALL wrk_dealloc( jkmax+4,jkmax+4, zrl01, kistart = 0, kjstart = 0 ) 
     148      CALL wrk_dealloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 ) 
     149      CALL wrk_dealloc( jpij, nlay_i+1, zqh_cum1, zh_cum1, kjstart = 0 ) 
     150      CALL wrk_dealloc( jpij, zhnew ) 
    692151      ! 
    693152   END SUBROUTINE lim_thd_ent 
Note: See TracChangeset for help on using the changeset viewer.