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 6625 for branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 – NEMO

Ignore:
Timestamp:
2016-05-26T11:08:07+02:00 (8 years ago)
Author:
kingr
Message:

Rolled back to r6613

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r6617 r6625  
    7474 
    7575      REAL(wp) ::   ztmelts             ! local scalar 
    76       REAL(wp) ::   zdum        
     76      REAL(wp) ::   zfdum        
    7777      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    7878      REAL(wp) ::   zs_snic      ! snow-ice salinity 
     
    9595      REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2) 
    9696      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2) 
    97       REAL(wp), POINTER, DIMENSION(:) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2) 
    9897 
    9998      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt  
     
    106105 
    107106      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2) 
     107      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2) 
     108      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3) 
    108109      REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing 
    109110 
     
    121122      END SELECT 
    122123 
    123       CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema ) 
    124       CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i ) 
     124      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 
     125      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
    125126      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 
    126127      CALL wrk_alloc( jpij, nlay_i, icount ) 
    127128        
    128       dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp ; dh_i_sub(:) = 0._wp 
     129      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp 
    129130      dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp    
    130131 
    131132      zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp 
    132       zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp ; zevap_rema(:) = 0._wp ; 
     133      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp 
    133134      zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp 
     135      zqh_s    (:) = 0._wp ; zq_s     (:) = 0._wp      
    134136 
    135137      zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp        
     
    157159      ! 
    158160      DO ji = kideb, kiut 
    159          zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
     161         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
    160162         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
    161163 
    162          zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
     164         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
    163165         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    164166      END DO 
     
    185187      !  2) Computing layer thicknesses and enthalpies.            ! 
    186188      !------------------------------------------------------------! 
     189      ! 
     190      DO jk = 1, nlay_s 
     191         DO ji = kideb, kiut 
     192            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 
     193         END DO 
     194      END DO 
    187195      ! 
    188196      DO jk = 1, nlay_i 
     
    267275      END DO 
    268276 
    269       !------------------------------ 
    270       ! 3.2 Sublimation (part1: snow)  
    271       !------------------------------ 
     277      !---------------------- 
     278      ! 3.2 Snow sublimation  
     279      !---------------------- 
    272280      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 
    273281      ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 
     282      ! clem comment: ice should also sublimate 
    274283      zdeltah(:,:) = 0._wp 
    275       DO ji = kideb, kiut 
    276          zdh_s_sub(ji)  = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
    277          ! remaining evap in kg.m-2 (used for ice melting later on) 
    278          zevap_rema(ji)  = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn 
    279          ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow) 
     284      ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice 
     285      ! forced  mode: snow thickness change due to sublimation 
     286      DO ji = kideb, kiut 
     287         zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
     288         ! Heat flux by sublimation [W.m-2], < 0 
     289         !      sublimate first snow that had fallen, then pre-existing snow 
    280290         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
    281291         hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  & 
     
    299309      !------------------------------------------- 
    300310      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 
     311      zq_s(:) = 0._wp  
    301312      DO jk = 1, nlay_s 
    302313         DO ji = kideb,kiut 
    303             rswitch       = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) ) 
    304             q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           & 
    305               &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  & 
    306               &              ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
     314            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  ) 
     315            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          & 
     316              &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  & 
     317              &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
     318            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk) 
    307319         END DO 
    308320      END DO 
     
    358370               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    359371                
    360                ! Contribution to salt flux >0 (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     372               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
    361373               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    362374                
     
    371383                
    372384            END IF 
    373             ! ---------------------- 
    374             ! Sublimation part2: ice 
    375             ! ---------------------- 
    376             zdum      = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic ) 
    377             zdeltah(ji,jk) = zdeltah(ji,jk) + zdum 
    378             dh_i_sub(ji)  = dh_i_sub(ji) + zdum 
    379             ! Salt flux > 0 (clem2016: flux is sent to the ocean for simplicity but salt should remain in the ice except if all ice is melted. 
    380             !                          It must be corrected at some point) 
    381             sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * sm_i_1d(ji) * r1_rdtice 
    382             ! Heat flux [W.m-2], < 0 
    383             hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice 
    384             ! Mass flux > 0 
    385             wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice 
    386             ! update remaining mass flux 
    387             zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic 
    388              
    389385            ! record which layers have disappeared (for bottom melting)  
    390386            !    => icount=0 : no layer has vanished 
     
    393389            icount(ji,jk) = NINT( rswitch ) 
    394390            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
    395                          
     391 
    396392            ! update heat content (J.m-2) and layer thickness 
    397393            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 
     
    401397      ! update ice thickness 
    402398      DO ji = kideb, kiut 
    403          ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) ) 
    404       END DO 
    405  
    406       ! remaining "potential" evap is sent to ocean 
    407       DO ji = kideb, kiut 
    408          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    409          wfx_err_sub(ii,ij) = wfx_err_sub(ii,ij) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1) 
     399         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) ) 
    410400      END DO 
    411401 
     
    696686      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 
    697687       
    698       CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema ) 
    699       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i ) 
     688      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 
     689      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
    700690      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 
    701691      CALL wrk_dealloc( jpij, nlay_i, icount ) 
Note: See TracChangeset for help on using the changeset viewer.