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 5123 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 – NEMO

Ignore:
Timestamp:
2015-03-04T17:06:03+01:00 (9 years ago)
Author:
clem
Message:

major LIM3 cleaning + monocat capabilities + NEMO namelist-consistency; sette to follow

File:
1 edited

Legend:

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

    r4990 r5123  
    2020   USE sbc_oce        ! Surface boundary condition: ocean fields 
    2121   USE ice            ! LIM variables 
    22    USE par_ice        ! LIM parameters 
    2322   USE thd_ice        ! LIM thermodynamics 
    2423   USE in_out_manager ! I/O manager 
     
    7069 
    7170      REAL(wp) ::   ztmelts             ! local scalar 
    72       REAL(wp) ::   zdh, zfdum  ! 
     71      REAL(wp) ::   zfdum        
    7372      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    7473      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
     
    9190      REAL(wp), POINTER, DIMENSION(:) ::   zq_su       ! heat for surface ablation                   (J.m-2) 
    9291      REAL(wp), POINTER, DIMENSION(:) ::   zq_bo       ! heat for bottom ablation                    (J.m-2) 
    93       REAL(wp), POINTER, DIMENSION(:) ::   zq_1cat     ! corrected heat in case 1-cat and hmelt>15cm (J.m-2) 
    9492      REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2) 
    95       REAL(wp), POINTER, DIMENSION(:) ::   zf_tt     ! Heat budget to determine melting or freezing(W.m-2) 
     93      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2) 
    9694      INTEGER , POINTER, DIMENSION(:) ::   icount      ! number of layers vanished by melting  
    9795 
     
    107105      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3) 
    108106 
    109       ! mass and salt flux (clem) 
    110       REAL(wp) :: zdvres, zswitch_sal 
     107      REAL(wp) :: zswitch_sal 
    111108 
    112109      ! Heat conservation  
     
    115112      !!------------------------------------------------------------------ 
    116113 
    117       ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values) 
    118       SELECT CASE( num_sal )                       ! varying salinity or not 
     114      ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values) 
     115      SELECT CASE( nn_icesal )                       ! varying salinity or not 
    119116         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile 
    120117         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile 
    121118      END SELECT 
    122119 
    123       CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 
     120      CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
    124121      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
    125122      CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) 
     
    130127  
    131128      zqprec (:) = 0._wp ; zq_su  (:) = 0._wp ; zq_bo  (:) = 0._wp ; zf_tt  (:) = 0._wp 
    132       zq_1cat(:) = 0._wp ; zq_rema(:) = 0._wp 
     129      zq_rema(:) = 0._wp 
    133130 
    134131      zh_s     (:) = 0._wp        
     
    148145      DO jk = 1, nlay_i 
    149146         DO ji = kideb, kiut 
    150             h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     147            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 
    151148            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 
    152149         ENDDO 
     
    159156      DO ji = kideb, kiut 
    160157         rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
    161          ztmelts       = rswitch * rtt + ( 1._wp - rswitch ) * rtt 
     158         ztmelts       = rswitch * rt0 + ( 1._wp - rswitch ) * rt0 
    162159 
    163160         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
     
    174171      !------------------------------------------------------------------------------! 
    175172      DO ji = kideb, kiut 
    176          IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting 
     173         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 
    177174            ! Contribution to heat flux to the ocean [W.m-2], < 0   
    178175            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 
     
    182179            ht_s_1d(ji)   = 0._wp 
    183180            q_s_1d (ji,1) = 0._wp 
    184             t_s_1d (ji,1) = rtt 
     181            t_s_1d (ji,1) = rt0 
    185182         END IF 
    186183      END DO 
     
    191188      ! 
    192189      DO ji = kideb, kiut      
    193          zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 
     190         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 
    194191      END DO 
    195192      ! 
     
    202199      DO jk = 1, nlay_i 
    203200         DO ji = kideb, kiut 
    204             zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     201            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 
    205202            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 
    206203         END DO 
     
    230227         !----------- 
    231228         ! thickness change 
    232          zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji)  
    233          zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 
     229         zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji)  
     230         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn 
    234231         ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 
    235          zqprec   (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )    
     232         zqprec   (ji) = rhosn * ( cpic * ( rt0 - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )    
    236233         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 
    237234         ! heat flux from snow precip (>0, W.m-2) 
     
    258255         zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) )       
    259256         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 
    260          zh_s  (ji) = ht_s_1d(ji) / REAL( nlay_s ) 
     257         zh_s  (ji) = ht_s_1d(ji) * r1_nlay_s 
    261258 
    262259         ENDIF 
     
    279276 
    280277            ! updates available heat + thickness 
    281             zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
     278            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
    282279            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 
    283280 
     
    314311      DO ji = kideb, kiut 
    315312         dh_s_tot(ji)   = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 
    316          zh_s(ji)       = ht_s_1d(ji) / REAL( nlay_s ) 
     313         zh_s(ji)       = ht_s_1d(ji) * r1_nlay_s 
    317314      END DO ! ji 
    318315 
     
    327324            q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) *             & 
    328325              &            ( (   MAX( 0._wp, dh_s_tot(ji) )              ) * zqprec(ji) +  & 
    329               &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) ) 
     326              &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
    330327            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk) 
    331328         END DO 
     
    338335      DO jk = 1, nlay_i 
    339336         DO ji = kideb, kiut  
    340             zEi            = - q_i_1d(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0] 
    341  
    342             ztmelts        = - tmut * s_i_1d(ji,jk) + rtt           ! Melting point of layer k [K] 
     337            zEi            = - q_i_1d(ji,jk) * r1_rhoic             ! Specific enthalpy of layer k [J/kg, <0] 
     338 
     339            ztmelts        = - tmut * s_i_1d(ji,jk) + rt0           ! Melting point of layer k [K] 
    343340 
    344341            zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
     
    348345            zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
    349346 
    350             zdeltah(ji,jk) = - zfmdt / rhoic                       ! Melt of layer jk [m, <0] 
     347            zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
    351348 
    352349            zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     
    408405      ! -> need for an iterative procedure, which converges quickly 
    409406 
    410       IF ( num_sal == 2 ) THEN 
     407      IF ( nn_icesal == 2 ) THEN 
    411408         num_iter_max = 5 
    412409      ELSE 
     
    414411      ENDIF 
    415412 
    416       !clem debug. Just to be sure that enthalpy at nlay_i+1 is null 
     413      ! Just to be sure that enthalpy at nlay_i+1 is null 
    417414      DO ji = kideb, kiut 
    418415         q_i_1d(ji,nlay_i+1) = 0._wp 
     
    440437                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji)  
    441438               ! New ice growth 
    442                ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     439               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) 
    443440 
    444441               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    445442                
    446443               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
    447                   &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   & 
    448                   &               + rcp  * ( ztmelts-rtt )           
     444                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & 
     445                  &               + rcp  * ( ztmelts-rt0 )           
    449446 
    450447               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     
    467464            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0) 
    468465 
    469             ztmelts        = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     466            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) 
    470467             
    471468            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    472469             
    473470            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
    474                &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   & 
    475                &               + rcp  * ( ztmelts-rtt )           
     471               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & 
     472               &               + rcp  * ( ztmelts-rt0 )           
    476473             
    477474            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     
    503500      DO jk = nlay_i, 1, -1 
    504501         DO ji = kideb, kiut 
    505             IF(  zf_tt(ji)  >=  0._wp  .AND. jk > icount(ji) ) THEN   ! do not calculate where layer has already disappeared from surface melting  
    506  
    507                ztmelts = - tmut * s_i_1d(ji,jk) + rtt  ! Melting point of layer jk (K) 
     502            IF(  zf_tt(ji)  >=  0._wp  .AND. jk > icount(ji) ) THEN   ! do not calculate where layer has already disappeared by surface melting  
     503 
     504               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K) 
    508505 
    509506               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
    510507 
    511                   zEi               = - q_i_1d(ji,jk) / rhoic        ! Specific enthalpy of melting ice (J/kg, <0) 
    512  
    513                   !!zEw               = rcp * ( t_i_1d(ji,jk) - rtt )  ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 
     508                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
     509 
     510                  !!zEw               = rcp * ( t_i_1d(ji,jk) - rt0 )  ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 
    514511 
    515512                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
     
    538535               ELSE                               !!! Basal melting 
    539536 
    540                   zEi               = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
    541  
    542                   zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) 
    543  
    544                   zdE               = zEi - zEw              ! Specific enthalpy difference   (J/kg, <0) 
    545  
    546                   zfmdt             = - zq_bo(ji) / zdE  ! Mass flux x time step (kg/m2, >0) 
    547  
    548                   zdeltah(ji,jk)    = - zfmdt / rhoic        ! Gross thickness change 
     537                  zEi               = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
     538 
     539                  zEw               = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0) 
     540 
     541                  zdE               = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
     542 
     543                  zfmdt             = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0) 
     544 
     545                  zdeltah(ji,jk)    = - zfmdt * r1_rhoic         ! Gross thickness change 
    549546 
    550547                  zdeltah(ji,jk)    = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
     
    576573            
    577574            ENDIF 
    578          END DO ! ji 
    579       END DO ! jk 
    580  
    581       !------------------------------------------------------------------------------! 
    582       ! Excessive ablation in a 1-category model 
    583       !     in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) 
    584       !------------------------------------------------------------------------------! 
    585       ! ??? keep ??? 
    586       ! clem bug: I think this should be included above, so we would not have to  
    587       !           track heat/salt/mass fluxes backwards 
    588 !      IF( jpl == 1 ) THEN 
    589 !         DO ji = kideb, kiut 
    590 !            IF(  zf_tt(ji)  >=  0._wp  ) THEN 
    591 !               zdh            = MAX( hmelt , dh_i_bott(ji) ) 
    592 !               zdvres         = zdh - dh_i_bott(ji) ! >=0 
    593 !               dh_i_bott(ji)  = zdh 
    594 ! 
    595 !               ! excessive energy is sent to lateral ablation 
    596 !               rswitch = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) ) 
    597 !               zq_1cat(ji) =  rswitch * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0 
    598 ! 
    599 !               ! correct salt and mass fluxes 
    600 !               sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation 
    601 !               wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdvres * r1_rdtice 
    602 !            ENDIF 
    603 !         END DO 
    604 !      ENDIF 
     575         END DO 
     576      END DO 
    605577 
    606578      !------------------------------------------- 
     
    635607         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    636608         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 
    637          hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
     609         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
    638610 
    639611         IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 
     
    655627         ! Salinity of snow ice 
    656628         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    657          zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
     629         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
    658630 
    659631         ! entrapment during snow ice formation 
    660632         ! new salinity difference stored (to be used in limthd_ent.F90) 
    661          IF (  num_sal == 2  ) THEN 
     633         IF (  nn_icesal == 2  ) THEN 
    662634            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 
    663635            ! salinity dif due to snow-ice formation 
     
    703675      DO ji = kideb, kiut 
    704676         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )  
    705          t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rtt 
     677         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 
    706678      END DO  ! ji 
    707679 
     
    712684            q_s_1d(ji,jk) = ( 1.0 - rswitch ) * q_s_1d(ji,jk) 
    713685            ! recalculate t_s_1d from q_s_1d 
    714             t_s_1d(ji,jk) = rtt + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 
     686            t_s_1d(ji,jk) = rt0 + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 
    715687         END DO 
    716688      END DO 
    717  
    718       CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 
     689       
     690      CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
    719691      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
    720692      CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) 
Note: See TracChangeset for help on using the changeset viewer.