Ignore:
Timestamp:
2013-11-07T11:01:27+01:00 (8 years ago)
Author:
cetlod
Message:

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r3808 r4161  
    66   !! History :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D 
    77   !!                 ! 2005-06 (M. Vancoppenolle) 3D version  
    8    !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif & rdmicif 
     8   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw & rdm_ice 
    99   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
    1010   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes  
     
    3939 
    4040   !!---------------------------------------------------------------------- 
    41    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     41   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    4242   !! $Id$ 
    4343   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    7373      !!  
    7474      INTEGER  ::   ji , jk        ! dummy loop indices 
    75       INTEGER  ::   ii, ij       ! 2D corresponding indices to ji 
     75      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji 
    7676      INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow 
    7777      INTEGER  ::   isnowic        ! snow ice formation not 
     
    8484      REAL(wp) ::   zdhnm, zhnnew, zhisn, zihic, zzc       ! 
    8585      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    86       REAL(wp) ::   zds          ! increment of bottom ice salinity 
    8786      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
    8887      REAL(wp) ::   zsm_snowice  ! snow-ice salinity 
     
    107106      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre   ! snow precipitation  
    108107      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub   ! snow sublimation 
    109       REAL(wp), POINTER, DIMENSION(:) ::   zsfx_melt   ! salt flux due to ice melt 
    110108 
    111109      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah 
     
    120118      REAL(wp), POINTER, DIMENSION(:,:) ::   zqt_i_lay   ! total ice heat content 
    121119 
     120      ! mass and salt flux (clem) 
     121      REAL(wp) :: zdvres, zdvsur, zdvbot 
     122      REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume... 
     123 
    122124      ! Heat conservation  
    123125      INTEGER  ::   num_iter_max, numce_dh 
     
    128130 
    129131      CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 
    130       CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zsfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
     132      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
    131133      CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i ) 
    132134      CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) 
    133135 
    134       zsfx_melt (:) = 0._wp 
     136      CALL wrk_alloc( jpij, zviold, zvsold ) ! clem 
     137       
    135138      ftotal_fin(:) = 0._wp 
    136139      zfdt_init (:) = 0._wp 
    137140      zfdt_final(:) = 0._wp 
    138141 
     142      dh_i_surf (:) = 0._wp 
     143      dh_i_bott (:) = 0._wp 
     144      dh_snowice(:) = 0._wp 
     145 
    139146      DO ji = kideb, kiut 
    140147         old_ht_i_b(ji) = ht_i_b(ji) 
    141148         old_ht_s_b(ji) = ht_s_b(ji) 
     149         zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem 
     150         zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem 
    142151      END DO 
    143152      ! 
     
    164173      ! 
    165174      DO ji = kideb, kiut     ! Layer thickness 
    166          zh_i(ji) = ht_i_b(ji) / nlay_i 
    167          zh_s(ji) = ht_s_b(ji) / nlay_s 
     175         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 
     176         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 
    168177      END DO 
    169178      ! 
     
    171180      DO jk = 1, nlay_s 
    172181         DO ji = kideb, kiut 
    173             zqt_s(ji) =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 
     182            zqt_s(ji) =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) 
    174183         END DO 
    175184      END DO 
     
    178187      DO jk = 1, nlay_i 
    179188         DO ji = kideb, kiut 
    180             zzc = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 
     189            zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 
    181190            zqt_i(ji)        =  zqt_i(ji) + zzc 
    182191            zqt_i_lay(ji,jk) =              zzc 
     
    244253         zhn            =  1.0 - MAX(  zzero , SIGN( zone , - zhsnew )  ) 
    245254         ht_s_b(ji)     =  MAX( zzero , zhsnew ) 
     255         ! we recompute dh_s_tot (clem)  
     256         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji) 
    246257         ! Volume and mass variations of snow 
    247258         dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) 
    248259         dvsbq_1d  (ji) =  MIN( zzero, dvsbq_1d(ji) ) 
    249          rdm_snw_1d(ji) =  rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 
     260         !clem rdm_snw_1d(ji) =  rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 
    250261      END DO ! ji 
    251262 
     
    254265      !-------------------------- 
    255266      DO ji = kideb, kiut  
    256          dh_i_surf(ji) =  0._wp 
    257267         z_f_surf (ji) =  zqfont_su(ji) * r1_rdtice   ! heat conservation test 
    258268         zdq_i    (ji) =  0._wp 
     
    272282            zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    273283            ! 
    274             !                                                    ! contribution to ice-ocean salt flux  
    275             zsfx_melt(ji)  = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice  
     284            ! clem 
     285            sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
     286               &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 
    276287         END DO 
    277288      END DO 
     
    334345         DO ji = kideb,kiut 
    335346            q_s_b    (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 
    336             zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s            ! heat conservation 
     347            zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s )            ! heat conservation 
    337348         END DO 
    338349      END DO 
     
    375386               ! Basal growth rate = - F*dt / q 
    376387               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
     388               sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 
    377389            ENDIF 
    378390         END DO 
     
    416428                  zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
    417429                     &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  
    418                   zds         = zfracs * sss_m(ii,ij) - s_i_new(ji) 
     430                  zfracs = MIN( 0.5 , zfracs ) 
    419431                  s_i_new(ji) = zfracs * sss_m(ii,ij) 
    420432               ENDIF ! fc_bo_i 
     
    425437         DO ji = kideb, kiut 
    426438            IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN 
    427                ! New ice salinity must not exceed 15 psu 
     439               ! New ice salinity must not exceed 20 psu 
    428440               s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 
    429441               ! Metling point in K 
     
    437449               ! Salinity update 
    438450               ! entrapment during bottom growth 
    439                dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) )    & 
    440                   &            / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji) 
     451               sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 
    441452            ENDIF ! heat budget 
    442453         END DO 
     
    476487                  zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    477488               ENDIF 
    478                ! contribution to salt flux 
    479                zsfx_melt(ji) = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice  
     489               ! clem: contribution to salt flux 
     490               sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
     491                    &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 
    480492            ENDIF 
    481493         END DO ! ji 
     
    528540         ELSE                  ;   zdhbf =              dh_i_bott(ji)  
    529541         ENDIF 
     542         zdvres        = zdhbf - dh_i_bott(ji) 
     543         dh_i_bott(ji) = zdhbf 
     544         sfx_thd_1d(ji)  = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 
    530545         !                     ! excessive energy is sent to lateral ablation 
    531          fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 )   & 
    532             &                          * ( zdhbf - dh_i_bott(ji) ) * r1_rdtice 
    533          dh_i_bott(ji)  = zdhbf 
    534          !                     !since ice volume is only used for outputs, we keep it global for all categories 
    535          dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 
    536          !                     !new ice thickness 
    537          zhgnew   (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 
    538          !                     ! diagnostic ( bottom ice growth ) 
    539          ii = MOD( npb(ji) - 1, jpi ) + 1 
    540          ij = ( npb(ji) - 1 ) / jpi + 1 
    541          diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
    542          diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 
    543          diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
     546         fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres * r1_rdtice 
    544547      END DO 
    545548 
     
    552555         ! Adapt the remaining energy if too much ice melts 
    553556         !-------------------------------------------------- 
    554          zihgnew =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice 
    555          ! 0 if no more ice 
    556          zhgnew    (ji) =         zihgnew   * zhgnew(ji)      ! ice thickness is put to 0 
    557          ! remaining heat 
     557         zdvres     = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) ) 
     558         zdvsur     = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first 
     559         zdvbot     = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom 
     560         dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem 
     561         dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem 
     562 
     563         ! new ice thickness (clem) 
     564         zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 
     565         zihgnew    = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 
     566         zhgnew(ji) = zihgnew * zhgnew(ji)      ! ice thickness is put to 0 
     567  
     568         !                     !since ice volume is only used for outputs, we keep it global for all categories 
     569         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 
     570 
     571        ! remaining heat 
    558572         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) )  
    559573 
     
    569583         ht_s_b(ji)     =  MAX( zzero , zhnfi ) 
    570584         zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji) 
     585         ! we recompute dh_s_tot (clem) 
     586         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji) 
    571587 
    572588         ! Mass variations of ice and snow 
     
    579595         ! 
    580596         !                                              ! mass variation cumulated over category 
    581          rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow  
    582          rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice  
     597         !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow  
     598         !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice  
    583599 
    584600         ! Remaining heat to the ocean  
     
    586602         focea(ji)  = - zfdt_final(ji) * r1_rdtice         ! focea is in W.m-2 * dt 
    587603 
     604         ! residual salt flux (clem) 
     605         !-------------------------- 
     606         ! surface 
     607         sfx_thd_1d(ji)    = sfx_thd_1d(ji) - sm_i_b(ji)  * a_i_b(ji) * zdvsur * rhoic * r1_rdtice 
     608         ! bottom 
     609         IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting 
     610            sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji)  * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 
     611         ELSE                                                          ! growth 
     612            sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 
     613         ENDIF 
     614         ! 
     615         ! diagnostic ( bottom ice growth ) 
     616         ii = MOD( npb(ji) - 1, jpi ) + 1 
     617         ij = ( npb(ji) - 1 ) / jpi + 1 
     618         diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
     619         diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 
     620         diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
    588621      END DO 
    589622 
     
    591624 
    592625      !--------------------------- 
    593       ! Salt flux and heat fluxes                     
     626      ! heat fluxes                     
    594627      !--------------------------- 
    595628      DO ji = kideb, kiut 
    596629         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice 
    597          ! 
    598          ! Salt flux 
    599          sfx_thd_1d(ji) = sfx_thd_1d(ji) +        zihgnew  * zsfx_melt(ji)               & 
    600             &                            - (1.0 - zihgnew) * zfmass_i (ji) * sm_i_b(ji)  * r1_rdtice 
    601630         ! 
    602631         ! Heat flux 
     
    646675         dmgwi_1d  (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 
    647676 
    648          ! All snow is thrown in the ocean, and seawater is taken to replace the volume 
    649          rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic * ( 1. - rhosn / rhoic ) 
    650          rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn 
     677         !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic  
     678         !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn  
    651679 
    652680         !        Equivalent salt flux (1) Snow-ice formation component 
     
    658686         ELSE                      ;   zsm_snowice = sm_i_b(ji)    
    659687         ENDIF 
    660          sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
    661          ! 
    662688         ! entrapment during snow ice formation 
    663          i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 
    664          isnowic      = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji)      ) ) * i_ice_switch 
    665          IF(  num_sal == 2  )   & 
    666             dsm_i_si_1d(ji) = (  zsm_snowice * dh_snowice(ji)    & 
    667             &                  + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 )   & 
    668             &                  - sm_i_b(ji)  ) * isnowic      
     689         ! clem: new salinity difference stored (to be used in limthd_ent.F90) 
     690         IF (  num_sal == 2  ) THEN 
     691            i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) ) 
     692            ! salinity dif due to snow-ice formation 
     693            dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch      
     694            ! salinity dif due to bottom growth  
     695            IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN 
     696               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 
     697            ENDIF 
     698         ENDIF 
    669699 
    670700         !  Actualize new snow and ice thickness. 
     
    680710         diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 
    681711         ! 
     712         ! salt flux 
     713         sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
     714         !-------------------------------- 
     715         ! Update mass fluxes (clem) 
     716         !-------------------------------- 
     717         rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic  
     718         rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn  
     719 
    682720      END DO !ji 
    683721      ! 
    684722      CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 
    685       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zsfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
     723      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
    686724      CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) 
    687725      CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) 
     726      ! 
     727      CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem 
    688728      ! 
    689729   END SUBROUTINE lim_thd_dh 
Note: See TracChangeset for help on using the changeset viewer.