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 15374 – NEMO

Changeset 15374


Ignore:
Timestamp:
2021-10-14T19:49:16+02:00 (10 months ago)
Author:
clem
Message:

rearrange (slightly) ice thermo routine

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icethd.F90

    r15334 r15374  
    9292      ! 
    9393      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices 
    94       REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos 
    95       REAL(wp), PARAMETER :: zfric_umin = 0._wp       ! lower bound for the friction velocity (cice value=5.e-04) 
    96       REAL(wp), PARAMETER :: zch        = 0.0057_wp   ! heat transfer coefficient 
    97       REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric, zvel   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 
    98       ! 
    99       ! for collection thickness 
    100       INTEGER  ::   iter             !   -       - 
    101       REAL(wp) ::   zvfrx, zvgx, ztaux, zf                     !   -      - 
    102       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp   !   -      - 
    103       REAL(wp), PARAMETER ::   zcai = 1.4e-3_wp               ! ice-air drag (clem: should be dependent on coupling/forcing used) 
    104       REAL(wp), PARAMETER ::   zhicrit = 0.04_wp                                             ! frazil ice thickness 
    105       REAL(wp), PARAMETER ::   zsqcd   = 1.0_wp / SQRT( 1.3_wp * zcai )                      ! 1/SQRT(airdensity*drag) 
    106       REAL(wp), PARAMETER ::   zgamafr = 0.03_wp 
    10794      !!------------------------------------------------------------------- 
    10895      ! controls 
     
    123110      ENDIF 
    124111 
    125       !---------------------------------------------! 
     112      ! --- calculate (some) heat fluxes and frazil ice --- ! 
     113      ! out => sensible heat (qsb_ice_bot) & heat budget in the leads (fhld & qlead) 
     114      ! out => ice collection thickness (ht_i_new) and fraction of frazil (fraz_frac) 
     115      CALL thd_prep 
     116       
     117      !-------------------------------------------------------------------------------------------! 
     118      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 
     119      !-------------------------------------------------------------------------------------------! 
     120      DO jl = 1, jpl 
     121 
     122         ! select ice covered grid points 
     123         npti = 0 ; nptidx(:) = 0 
     124         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     125            IF ( a_i(ji,jj,jl) > epsi10 ) THEN 
     126               npti         = npti  + 1 
     127               nptidx(npti) = (jj - 1) * jpi + ji 
     128            ENDIF 
     129         END_2D 
     130 
     131         IF( npti > 0 ) THEN  ! If there is no ice, do nothing. 
     132            ! 
     133                              CALL ice_thd_1d2d( jl, 1 )            ! --- Move to 1D arrays --- ! 
     134            !                                                       ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- ! 
     135            ! 
     136            s_i_new   (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp   ! --- some init --- !  (important to have them here) 
     137            dh_i_sum  (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm  (1:npti) = 0._wp 
     138            dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 
     139            dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 
     140            ! 
     141                              CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- ! 
     142            ! 
     143            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- ! 
     144                              CALL ice_thd_dh                           ! Ice-Snow thickness 
     145                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping 
     146            ENDIF 
     147                              CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- ! 
     148            ! 
     149                              CALL ice_thd_temp                     ! --- Temperature update --- ! 
     150            ! 
     151            IF( ln_icedH .AND. ln_virtual_itd ) & 
     152               &              CALL ice_thd_mono                     ! --- Extra lateral melting if virtual_itd --- ! 
     153            ! 
     154            IF( ln_icedA )    CALL ice_thd_da                       ! --- Lateral melting --- ! 
     155            ! 
     156                              CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- ! 
     157            !                                                       ! --- & Move to 2D arrays --- ! 
     158         ENDIF 
     159         ! 
     160      END DO 
     161      ! 
     162      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     163      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
     164      ! 
     165      IF ( ln_pnd .AND. ln_icedH ) & 
     166         &                    CALL ice_thd_pnd                      ! --- Melt ponds --- ! 
     167      ! 
     168      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
     169      ! 
     170      IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- ! 
     171      ! 
     172                              CALL ice_cor( kt , 2 )                ! --- Corrections --- ! 
     173      ! 
     174      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice              ! --- Ice natural aging incrementation 
     175      ! 
     176      DO_2D( 0, 0, 0, 0 )                                           ! --- Ice velocity corrections 
     177         IF( at_i(ji,jj) == 0._wp ) THEN   ! if ice has melted 
     178            IF( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side 
     179            IF( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp   ! left side 
     180            IF( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp   ! upper side 
     181            IF( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp   ! bottom side 
     182         ENDIF 
     183      END_2D 
     184      CALL lbc_lnk( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 
     185      ! 
     186      ! convergence tests 
     187      IF( ln_zdf_chkcvg ) THEN 
     188         CALL iom_put( 'tice_cvgerr', ztice_cvgerr ) ; DEALLOCATE( ztice_cvgerr ) 
     189         CALL iom_put( 'tice_cvgstp', ztice_cvgstp ) ; DEALLOCATE( ztice_cvgstp ) 
     190      ENDIF 
     191      ! 
     192      ! controls 
     193      IF( ln_icectl )   CALL ice_prt    (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints 
     194      IF( sn_cfctl%l_prtctl )   & 
     195        &               CALL ice_prt3D  ('icethd')                                        ! prints 
     196      IF( ln_timing )   CALL timing_stop('icethd')                                        ! timing 
     197      ! 
     198   END SUBROUTINE ice_thd 
     199 
     200 
     201   SUBROUTINE ice_thd_temp 
     202      !!----------------------------------------------------------------------- 
     203      !!                   ***  ROUTINE ice_thd_temp *** 
     204      !! 
     205      !! ** Purpose :   Computes sea ice temperature (Kelvin) from enthalpy 
     206      !! 
     207      !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
     208      !!------------------------------------------------------------------- 
     209      INTEGER  ::   ji, jk   ! dummy loop indices 
     210      REAL(wp) ::   ztmelts, zbbb, zccc  ! local scalar 
     211      !!------------------------------------------------------------------- 
     212      ! Recover ice temperature 
     213      DO jk = 1, nlay_i 
     214         DO ji = 1, npti 
     215            ztmelts       = -rTmlt * sz_i_1d(ji,jk) 
     216            ! Conversion q(S,T) -> T (second order equation) 
     217            zbbb          = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus 
     218            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) ) 
     219            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi 
     220 
     221            ! mask temperature 
     222            rswitch       = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
     223            t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0 
     224         END DO 
     225      END DO 
     226      ! 
     227   END SUBROUTINE ice_thd_temp 
     228 
     229 
     230   SUBROUTINE ice_thd_mono 
     231      !!----------------------------------------------------------------------- 
     232      !!                   ***  ROUTINE ice_thd_mono *** 
     233      !! 
     234      !! ** Purpose :   Lateral melting in case virtual_itd 
     235      !!                          ( dA = A/2h dh ) 
     236      !!----------------------------------------------------------------------- 
     237      INTEGER  ::   ji                 ! dummy loop indices 
     238      REAL(wp) ::   zhi_bef            ! ice thickness before thermo 
     239      REAL(wp) ::   zdh_mel, zda_mel   ! net melting 
     240      REAL(wp) ::   zvi, zvs           ! ice/snow volumes 
     241      !!----------------------------------------------------------------------- 
     242      ! 
     243      DO ji = 1, npti 
     244         zdh_mel = MIN( 0._wp, dh_i_itm(ji) + dh_i_sum(ji) + dh_i_bom(ji) + dh_snowice(ji) + dh_i_sub(ji) ) 
     245         IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN 
     246            zvi          = a_i_1d(ji) * h_i_1d(ji) 
     247            zvs          = a_i_1d(ji) * h_s_1d(ji) 
     248            ! lateral melting = concentration change 
     249            zhi_bef     = h_i_1d(ji) - zdh_mel 
     250            rswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) 
     251            zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 
     252            a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel ) 
     253            ! adjust thickness 
     254            h_i_1d(ji) = zvi / a_i_1d(ji) 
     255            h_s_1d(ji) = zvs / a_i_1d(ji) 
     256            ! retrieve total concentration 
     257            at_i_1d(ji) = a_i_1d(ji) 
     258         END IF 
     259      END DO 
     260      ! 
     261   END SUBROUTINE ice_thd_mono 
     262 
     263   SUBROUTINE thd_prep 
     264      !!----------------------------------------------------------------------- 
     265      !!                   ***  ROUTINE thd_prep *** 
     266      !! 
     267      !! ** Purpose :   prepare necessary fields for thermo calculations 
     268      !! 
     269      !! For the fluxes 
     270      !! ** Inputs  :   u_ice, v_ice, ssu_m, ssv_m, utau, vtau 
     271      !!                frq_m, qsr_oce, qns_oce, qemp_oce, e3t_m, sst_m 
     272      !! ** Outputs :   qsb_ice_bot, fhld, qlead 
     273      !! 
     274      !! For the collection thickness (frazil) 
     275      !! ** Inputs  :   u_ice, v_ice, utau_ice, vtau_ice 
     276      !! ** Ouputs  :   ht_i_new, fraz_frac 
     277      !!----------------------------------------------------------------------- 
     278      INTEGER  ::   ji, jj             ! dummy loop indices 
     279      REAL(wp) ::   zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos, zu_io, zv_io, zu_iom1, zv_iom1 
     280      REAL(wp), PARAMETER ::   zfric_umin = 0._wp       ! lower bound for the friction velocity (cice value=5.e-04) 
     281      REAL(wp), PARAMETER ::   zch        = 0.0057_wp   ! heat transfer coefficient 
     282      REAL(wp), DIMENSION(jpi,jpj) ::  zfric, zvel      ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 
     283      ! 
     284      ! for frazil ice 
     285      INTEGER  ::   iter 
     286      REAL(wp) ::   zvfrx, zvgx, ztaux, zf, ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp 
     287      REAL(wp), PARAMETER ::   zcai    = 1.4e-3_wp                       ! ice-air drag (clem: should be dependent on coupling/forcing used) 
     288      REAL(wp), PARAMETER ::   zhicrit = 0.04_wp                         ! frazil ice thickness 
     289      REAL(wp), PARAMETER ::   zsqcd   = 1.0_wp / SQRT( 1.3_wp * zcai )  ! 1/SQRT(airdensity*drag) 
     290      REAL(wp), PARAMETER ::   zgamafr = 0.03_wp 
     291      !!----------------------------------------------------------------------- 
     292      ! 
    126293      ! computation of friction velocity at T points 
    127       !---------------------------------------------! 
    128294      IF( ln_icedyn ) THEN 
    129          zu_io(:,:) = u_ice(:,:) - ssu_m(:,:) 
    130          zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 
    131295         DO_2D( 0, 0, 0, 0 ) 
    132             zfric(ji,jj) = rn_cio * ( 0.5_wp *  & 
    133                &                    (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
    134                &                     + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 
    135             zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + & 
    136                &                         ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 
     296            zu_io   = u_ice(ji  ,jj  ) - ssu_m(ji  ,jj  ) 
     297            zu_iom1 = u_ice(ji-1,jj  ) - ssu_m(ji-1,jj  ) 
     298            zv_io   = v_ice(ji  ,jj  ) - ssv_m(ji  ,jj  ) 
     299            zv_iom1 = v_ice(ji  ,jj-1) - ssv_m(ji  ,jj-1) 
     300            ! 
     301            zfric(ji,jj) = rn_cio * ( 0.5_wp * ( zu_io*zu_io + zu_iom1*zu_iom1 + zv_io*zv_io + zv_iom1*zv_iom1 ) ) * tmask(ji,jj,1) 
     302            zvel (ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj  ) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) ) + & 
     303               &                          ( v_ice(ji  ,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) ) ) 
    137304         END_2D 
    138305      ELSE      !  if no ice dynamics => transfer directly the atmospheric stress to the ocean 
     
    227394      ENDIF 
    228395 
    229  
    230396      !---------------------------------------------------------! 
    231397      ! Collection thickness of ice formed in leads and polynyas 
     
    251417            IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 
    252418               ! -- Wind stress -- ! 
    253                ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   & 
    254                   &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp 
    255                ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   & 
    256                   &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp 
     419               ztaux = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1) + utau_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 
     420               ztauy = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1) + vtau_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 
    257421               ! Square root of wind stress 
    258                ztenagm       = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
     422               ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
    259423 
    260424               ! -- Frazil ice velocity -- ! 
     
    268432 
    269433               ! -- Relative frazil/pack ice velocity -- ! 
    270                rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
    271                zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   & 
    272                   &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15_wp * 0.15_wp ) * rswitch 
    273  
    274                !!clem 
     434               rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
     435               zvrel2  = MAX( (zvfrx - zvgx)*(zvfrx - zvgx) + (zvfry - zvgy)*(zvfry - zvgy), 0.15_wp*0.15_wp ) * rswitch 
     436 
     437               ! -- fraction of frazil ice -- ! 
    275438               fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz 
    276                !!clem 
    277439                
    278440               ! -- new ice thickness (iterative loop) -- ! 
     
    301463 
    302464      ENDIF 
    303  
    304       !-------------------------------------------------------------------------------------------! 
    305       ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 
    306       !-------------------------------------------------------------------------------------------! 
    307       DO jl = 1, jpl 
    308  
    309          ! select ice covered grid points 
    310          npti = 0 ; nptidx(:) = 0 
    311          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    312             IF ( a_i(ji,jj,jl) > epsi10 ) THEN 
    313                npti         = npti  + 1 
    314                nptidx(npti) = (jj - 1) * jpi + ji 
    315             ENDIF 
    316          END_2D 
    317  
    318          IF( npti > 0 ) THEN  ! If there is no ice, do nothing. 
    319             ! 
    320                               CALL ice_thd_1d2d( jl, 1 )            ! --- Move to 1D arrays --- ! 
    321             !                                                       ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- ! 
    322             ! 
    323             s_i_new   (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp   ! --- some init --- !  (important to have them here) 
    324             dh_i_sum  (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm  (1:npti) = 0._wp 
    325             dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 
    326             dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 
    327             ! 
    328                               CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- ! 
    329             ! 
    330             IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- ! 
    331                               CALL ice_thd_dh                           ! Ice-Snow thickness 
    332                               CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping 
    333             ENDIF 
    334                               CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- ! 
    335             ! 
    336                               CALL ice_thd_temp                     ! --- Temperature update --- ! 
    337             ! 
    338             IF( ln_icedH .AND. ln_virtual_itd ) & 
    339                &              CALL ice_thd_mono                     ! --- Extra lateral melting if virtual_itd --- ! 
    340             ! 
    341             IF( ln_icedA )    CALL ice_thd_da                       ! --- Lateral melting --- ! 
    342             ! 
    343                               CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- ! 
    344             !                                                       ! --- & Move to 2D arrays --- ! 
    345          ENDIF 
    346          ! 
    347       END DO 
    348       ! 
    349       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    350       IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    351       ! 
    352       IF ( ln_pnd .AND. ln_icedH ) & 
    353          &                    CALL ice_thd_pnd                      ! --- Melt ponds --- ! 
    354       ! 
    355       IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
    356       ! 
    357       IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- ! 
    358       ! 
    359                               CALL ice_cor( kt , 2 )                ! --- Corrections --- ! 
    360       ! 
    361       oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice              ! --- Ice natural aging incrementation 
    362       ! 
    363       DO_2D( 0, 0, 0, 0 )                                           ! --- Ice velocity corrections 
    364          IF( at_i(ji,jj) == 0._wp ) THEN   ! if ice has melted 
    365             IF( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side 
    366             IF( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp   ! left side 
    367             IF( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp   ! upper side 
    368             IF( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp   ! bottom side 
    369          ENDIF 
    370       END_2D 
    371       CALL lbc_lnk( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 
    372       ! 
    373       ! convergence tests 
    374       IF( ln_zdf_chkcvg ) THEN 
    375          CALL iom_put( 'tice_cvgerr', ztice_cvgerr ) ; DEALLOCATE( ztice_cvgerr ) 
    376          CALL iom_put( 'tice_cvgstp', ztice_cvgstp ) ; DEALLOCATE( ztice_cvgstp ) 
    377       ENDIF 
    378       ! 
    379       ! controls 
    380       IF( ln_icectl )   CALL ice_prt    (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints 
    381       IF( sn_cfctl%l_prtctl )   & 
    382         &               CALL ice_prt3D  ('icethd')                                        ! prints 
    383       IF( ln_timing )   CALL timing_stop('icethd')                                        ! timing 
    384       ! 
    385    END SUBROUTINE ice_thd 
    386  
    387  
    388    SUBROUTINE ice_thd_temp 
    389       !!----------------------------------------------------------------------- 
    390       !!                   ***  ROUTINE ice_thd_temp *** 
    391       !! 
    392       !! ** Purpose :   Computes sea ice temperature (Kelvin) from enthalpy 
    393       !! 
    394       !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
    395       !!------------------------------------------------------------------- 
    396       INTEGER  ::   ji, jk   ! dummy loop indices 
    397       REAL(wp) ::   ztmelts, zbbb, zccc  ! local scalar 
    398       !!------------------------------------------------------------------- 
    399       ! Recover ice temperature 
    400       DO jk = 1, nlay_i 
    401          DO ji = 1, npti 
    402             ztmelts       = -rTmlt * sz_i_1d(ji,jk) 
    403             ! Conversion q(S,T) -> T (second order equation) 
    404             zbbb          = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus 
    405             zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) ) 
    406             t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi 
    407  
    408             ! mask temperature 
    409             rswitch       = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
    410             t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0 
    411          END DO 
    412       END DO 
    413       ! 
    414    END SUBROUTINE ice_thd_temp 
    415  
    416  
    417    SUBROUTINE ice_thd_mono 
    418       !!----------------------------------------------------------------------- 
    419       !!                   ***  ROUTINE ice_thd_mono *** 
    420       !! 
    421       !! ** Purpose :   Lateral melting in case virtual_itd 
    422       !!                          ( dA = A/2h dh ) 
    423       !!----------------------------------------------------------------------- 
    424       INTEGER  ::   ji                 ! dummy loop indices 
    425       REAL(wp) ::   zhi_bef            ! ice thickness before thermo 
    426       REAL(wp) ::   zdh_mel, zda_mel   ! net melting 
    427       REAL(wp) ::   zvi, zvs           ! ice/snow volumes 
    428       !!----------------------------------------------------------------------- 
    429       ! 
    430       DO ji = 1, npti 
    431          zdh_mel = MIN( 0._wp, dh_i_itm(ji) + dh_i_sum(ji) + dh_i_bom(ji) + dh_snowice(ji) + dh_i_sub(ji) ) 
    432          IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN 
    433             zvi          = a_i_1d(ji) * h_i_1d(ji) 
    434             zvs          = a_i_1d(ji) * h_s_1d(ji) 
    435             ! lateral melting = concentration change 
    436             zhi_bef     = h_i_1d(ji) - zdh_mel 
    437             rswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) 
    438             zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 
    439             a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel ) 
    440             ! adjust thickness 
    441             h_i_1d(ji) = zvi / a_i_1d(ji) 
    442             h_s_1d(ji) = zvs / a_i_1d(ji) 
    443             ! retrieve total concentration 
    444             at_i_1d(ji) = a_i_1d(ji) 
    445          END IF 
    446       END DO 
    447       ! 
    448    END SUBROUTINE ice_thd_mono 
    449  
     465       
     466   END SUBROUTINE thd_prep 
    450467 
    451468   SUBROUTINE ice_thd_1d2d( kl, kn ) 
Note: See TracChangeset for help on using the changeset viewer.