Changeset 2383


Ignore:
Timestamp:
2010-11-13T16:01:26+01:00 (10 years ago)
Author:
gm
Message:

v3.3beta: LIM-3 bug correction ticket #670

Location:
branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r2370 r2383  
    232232 
    233233            !  computing freshwater exchanges at the ice/ocean interface 
    234             zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj) )  &   !  evaporation over oceanic fraction 
    235                &   + tprecip(ji,jj) *         at_i(ji,jj)    &   !  total precipitation 
    236                ! old fashioned way                
    237                !              &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )  &   !  remov. snow precip over ice 
    238                &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   !  remov. snow precip over ice 
    239                &   - rdmsnif(ji,jj) * r1_rdtice                       &   !  freshwaterflux due to snow melting  
    240                ! new contribution from snow falling when ridging 
    241                &   + fmmec(ji,jj) 
     234            zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction 
     235               &   + tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean 
     236               &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice 
     237               &   - rdmsnif(ji,jj) * r1_rdtice                       &   ! freshwaterflux due to snow melting  
     238               &   + fmmec(ji,jj)                                         ! snow falling when ridging 
     239 
    242240 
    243241            !  computing salt exchanges at the ice/ocean interface 
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r2287 r2383  
    77   !!            2.0  ! 2002-07 (C. Ethe, G. Madec)  LIM-2 (F90 rewriting) 
    88   !!            3.0  ! 2005-11 (M. Vancoppenolle)  LIM-3 : Multi-layer thermodynamics + salinity variations 
    9    !!             -   ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec and lim_thd_con_dif 
     9   !!             -   ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec, lim_thd_con_dh and lim_thd_con_dif 
    1010   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif 
     11   !!            3.3  ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas) 
    1112   !!---------------------------------------------------------------------- 
    1213#if defined key_lim3 
     
    1415   !!   'key_lim3'                                      LIM3 sea-ice model 
    1516   !!---------------------------------------------------------------------- 
    16    !!   lim_thd      : thermodynamic of sea ice 
    17    !!   lim_thd_init : initialisation of sea-ice thermodynamic 
     17   !!   lim_thd        : thermodynamic of sea ice 
     18   !!   lim_thd_init   : initialisation of sea-ice thermodynamic 
    1819   !!---------------------------------------------------------------------- 
    1920   USE phycst          ! physical constants 
    2021   USE dom_oce         ! ocean space and time domain variables 
    21    USE ice             ! LIM sea-ice variables 
    22    USE par_ice 
     22   USE ice             ! LIM: sea-ice variables 
     23   USE par_ice         ! LIM: sea-ice parameters 
    2324   USE sbc_oce         ! Surface boundary condition: ocean fields 
    2425   USE sbc_ice         ! Surface boundary condition: ice fields 
    2526   USE thd_ice         ! LIM thermodynamic sea-ice variables 
    2627   USE dom_ice         ! LIM sea-ice domain 
    27    USE domvvl 
    28    USE iceini 
    29    USE limthd_dif 
    30    USE limthd_dh 
    31    USE limthd_sal 
    32    USE limthd_ent 
    33    USE limtab 
    34    USE limvar 
     28   USE domvvl          ! domain: variable volume level 
     29   USE iceini          ! LIM: sea-ice initialization 
     30   USE limthd_dif      ! LIM: thermodynamics, vertical diffusion 
     31   USE limthd_dh       ! LIM: thermodynamics, ice and snow thickness variation 
     32   USE limthd_sal      ! LIM: thermodynamics, ice salinity 
     33   USE limthd_ent      ! LIM: thermodynamics, ice enthalpy redistribution 
     34   USE limtab          ! LIM: 1D <==> 2D transformation 
     35   USE limvar          ! LIM: ??? 
     36   USE lbclnk          ! lateral boundary condition - MPP links 
     37   USE lib_mpp         ! MPP library 
    3538   USE in_out_manager  ! I/O manager 
    3639   USE prtctl          ! Print control 
    37    USE lbclnk 
    38    USE lib_mpp 
    3940 
    4041   IMPLICIT NONE 
     
    4344   PUBLIC   lim_thd    ! called by lim_step 
    4445 
    45    REAL(wp) ::   epsi20 = 1e-20   ! constant values 
    46    REAL(wp) ::   epsi16 = 1e-16   ! 
    47    REAL(wp) ::   epsi06 = 1e-06   ! 
    48    REAL(wp) ::   epsi04 = 1e-04   ! 
    49    REAL(wp) ::   zzero  = 0.e0    ! 
    50    REAL(wp) ::   zone   = 1.e0    ! 
     46   REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values 
     47   REAL(wp) ::   epsi16 = 1e-16_wp   ! 
     48   REAL(wp) ::   epsi06 = 1e-06_wp   ! 
     49   REAL(wp) ::   epsi04 = 1e-04_wp   ! 
     50   REAL(wp) ::   zzero  = 0._wp      ! 
     51   REAL(wp) ::   zone   = 1._wp      ! 
    5152 
    5253   !! * Substitutions 
     
    5657   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
    5758   !! $Id$ 
    58    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    5960   !!---------------------------------------------------------------------- 
    60  
    6161CONTAINS 
    6262 
     
    8686      REAL(wp) ::   zfric_umax = 2e-02    ! upper bound for the friction velocity 
    8787      REAL(wp) ::   zinda, zindb, zthsnice, zfric_u    ! temporary scalar 
    88       REAL(wp) ::   zfnsol, zfontn, zfntlat, zpareff   !    -         - 
     88      REAL(wp) ::   zfntlat, zpareff   !    -         - 
    8989      REAL(wp) ::   zeps, zareamin, zcoef 
    9090      REAL(wp), DIMENSION(jpi,jpj) ::   zqlbsbq   ! link with lead energy budget qldif 
     
    104104      !-------------------- 
    105105      ! Change the units of heat content; from global units to J.m3 
    106  
    107106      DO jl = 1, jpl 
    108107         DO jk = 1, nlay_i 
     
    171170            pfrld(ji,jj)   = 1.0 - at_i(ji,jj) 
    172171            zinda          = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) ) 
    173  
     172            ! 
    174173            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
    175174            !           !  practically no "direct lateral ablation" 
     
    186185            qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 
    187186            !                        
    188  
    189             ! still need to be updated : fdtcn !!!! 
    190             !           !-- Lead heat budget (part 1, next one is in limthd_dh 
    191             !           !-- qldif -- (or qldif_1d in 1d routines) 
    192             zfontn         = sprecip(ji,jj) * lfus              ! energy of melting 
    193             zfnsol         = qns(ji,jj)                         ! total non solar flux 
    194             qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj)                               & 
    195                &                               + zfnsol + fdtcn(ji,jj) - zfontn      & 
    196                &                               + ( 1.0 - zindb ) * fsbbq(ji,jj)  )   & 
    197                &                               * ( 1.0 - at_i(ji,jj) ) * rdt_ice     
    198  
     187            !           !-- Lead heat budget, qldif (part 1, next one is in limthd_dh)  
     188            !           !   caution: exponent betas used as more snow can fallinto leads 
     189            qldif(ji,jj) =  tms(ji,jj) * rdt_ice  * (                             & 
     190               &   pfrld(ji,jj)        * (  qsr(ji,jj)                            &   ! solar heat 
     191               &                            + qns(ji,jj)                          &   ! non solar heat 
     192               &                            + fdtcn(ji,jj)                        &   ! turbulent ice-ocean heat 
     193               &                            + fsbbq(ji,jj) * ( 1.0 - zindb )  )   &   ! residual heat from previous step 
     194               & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus                    )   ! latent heat of sprecip melting 
     195            ! 
    199196            ! Positive heat budget is used for bottom ablation 
    200197            zfntlat        = 1.0 - MAX( zzero , SIGN( zone ,  - qldif(ji,jj) ) ) 
     
    203200            != 0 if ice and positive heat budget and 1 if one of those two is false 
    204201            zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( at_i(ji,jj) * rdt_ice , epsi16 ) 
    205  
     202            ! 
    206203            ! Heat budget of the lead, energy transferred from ice to ocean 
    207204            qldif  (ji,jj) = zpareff * qldif(ji,jj) 
    208205            qdtcn  (ji,jj) = zpareff * qdtcn(ji,jj) 
    209  
     206            ! 
    210207            ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 
    211208            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda ) 
    212  
     209            ! 
    213210            ! oceanic heat flux (limthd_dh) 
    214211            fbif   (ji,jj) = zindb * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 
     
    228225         ENDIF 
    229226 
    230          zareamin = 1.0e-10 
     227         zareamin = 1.e-10 
    231228         nbpb = 0 
    232229         DO jj = 1, jpj 
     
    364361            CALL tab_1d_2d( nbpb, rdvonif, npb, dvnbq_1d  (1:nbpb), jpi, jpj )  
    365362            CALL tab_1d_2d( nbpb, fseqv  , npb, fseqv_1d  (1:nbpb), jpi, jpj ) 
    366  
     363            ! 
    367364            IF( num_sal == 2 ) THEN 
    368365               CALL tab_1d_2d( nbpb, fsbri, npb, fsbri_1d(1:nbpb), jpi, jpj ) 
    369366               CALL tab_1d_2d( nbpb, fhbri, npb, fhbri_1d(1:nbpb), jpi, jpj ) 
    370367            ENDIF 
    371  
     368            ! 
    372369            !+++++ 
    373370            !temporary stuff for a dummy version 
     
    380377            CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj) 
    381378            !+++++ 
    382  
     379            ! 
    383380            IF( lk_mpp )   CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 
    384381         ENDIF 
     
    393390      ! 5.1) Ice heat content               
    394391      !------------------------ 
    395  
    396392      ! Enthalpies are global variables we have to readjust the units 
    397393      zcoef = 1.e0 / ( unit_fac * REAL(nlay_i) ) 
     
    406402      ! 5.2) Snow heat content               
    407403      !------------------------ 
    408  
    409404      ! Enthalpies are global variables we have to readjust the units 
    410405      zcoef = 1.e0 / ( unit_fac * REAL(nlay_s) ) 
     
    429424      IF( con_i )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 
    430425 
    431       IF(ln_ctl) THEN   ! Control print 
     426      IF(ln_ctl) THEN            ! Control print 
    432427         CALL prt_ctl_info(' ') 
    433428         CALL prt_ctl_info(' - Cell values : ') 
     
    459454            END DO 
    460455         END DO 
    461  
    462456      ENDIF 
    463  
     457      ! 
    464458   END SUBROUTINE lim_thd 
    465459 
     
    529523      END DO 
    530524      ! layer by layer 
    531       dq_i_layer(:,:)    = q_i_layer_fin(:,:) - q_i_layer_in(:,:) 
     525      dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:) 
    532526 
    533527      !---------------------------------------- 
    534528      ! Atmospheric heat flux, ice heat budget 
    535529      !---------------------------------------- 
    536  
    537       DO ji = kideb, kiut 
    538          zji = MOD( npb(ji) - 1, jpi ) + 1 
    539          zjj = ( npb(ji) - 1 ) / jpi + 1 
    540  
    541          fatm(ji,jl) = qnsr_ice_1d(ji) + (1.0-i0(ji))*qsr_ice_1d(ji) 
    542  
    543          sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji)*i0(ji) - fstroc(zji,zjj,jl) 
     530      DO ji = kideb, kiut 
     531         zji = MOD( npb(ji) - 1 , jpi ) + 1 
     532         zjj =    ( npb(ji) - 1 ) / jpi + 1 
     533         fatm     (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) 
     534         sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(zji,zjj,jl) 
    544535      END DO 
    545536 
     
    547538      ! Conservation error 
    548539      !-------------------- 
    549  
    550540      DO ji = kideb, kiut 
    551541         cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 
    552542      END DO 
    553543 
    554       numce = 0 
     544      numce  = 0 
    555545      meance = 0.0 
    556546      DO ji = kideb, kiut 
     
    559549            meance = meance + cons_error(ji,jl) 
    560550         ENDIF 
    561       ENDDO 
    562       IF (numce .GT. 0 ) meance = meance / numce 
     551      END DO 
     552      IF( numce .GT. 0 ) meance = meance / numce 
    563553 
    564554      WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err 
     
    570560      ! Surface error due to imbalance between Fatm and Fcsu 
    571561      !------------------------------------------------------- 
    572       numce  = 0.0 
     562      numce  = 0 
    573563      meance = 0.0 
    574564 
    575565      DO ji = kideb, kiut 
    576566         surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) ) 
    577          IF ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. & 
    578             max_surf_err ) ) THEN 
     567         IF( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) THEN 
    579568            numce = numce + 1  
    580569            meance = meance + surf_error(ji,jl) 
    581570         ENDIF 
    582571      ENDDO 
    583       IF (numce .GT. 0 ) meance = meance / numce 
     572      IF( numce .GT. 0 ) meance = meance / numce 
    584573 
    585574      WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err 
     
    595584      ! Write ice state in case of big errors 
    596585      !--------------------------------------- 
    597  
    598586      DO ji = kideb, kiut 
    599587         IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. & 
     
    601589            zji                 = MOD( npb(ji) - 1, jpi ) + 1 
    602590            zjj                 = ( npb(ji) - 1 ) / jpi + 1 
    603  
     591            ! 
    604592            WRITE(numout,*) ' alerte 1     ' 
    605593            WRITE(numout,*) ' Untolerated conservation / surface error after ' 
     
    617605            !        WRITE(numout,*) ' qt_i_fin  : ', qt_i_fin(ji,jl) 
    618606            !        WRITE(numout,*) ' qt_s_fin  : ', qt_s_fin(ji,jl) 
    619             !        WRITE(numout,*) ' qt        : ', qt_i_fin(ji,jl) + & 
    620             !                                         qt_s_fin(ji,jl) 
     607            !        WRITE(numout,*) ' qt        : ', qt_i_fin(ji,jl) + qt_s_fin(ji,jl) 
    621608            WRITE(numout,*) ' ht_i       : ', ht_i_b(ji) 
    622609            WRITE(numout,*) ' ht_s       : ', ht_s_b(ji) 
     
    645632            WRITE(numout,*) 
    646633            WRITE(numout,*) ' Layer by layer ... ' 
    647             WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - & 
    648                qt_s_in(ji,jl) )  &  
    649                / rdt_ice 
    650             WRITE(numout,*) ' dfc_snow  : ', fc_s(ji,1) -      & 
    651                fc_s(ji,0) 
     634            WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice 
     635            WRITE(numout,*) ' dfc_snow  : ', fc_s(ji,1) - fc_s(ji,0) 
    652636            DO jk = 1, nlay_i 
    653637               WRITE(numout,*) ' layer  : ', jk 
    654638               WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice   
    655639               WRITE(numout,*) ' radab  : ', radab(ji,jk) 
    656                WRITE(numout,*) ' dfc_i  : ', fc_i(ji,jk) -               & 
    657                   fc_i(ji,jk-1) 
    658                WRITE(numout,*) ' tot f  : ', fc_i(ji,jk) -               & 
    659                   fc_i(ji,jk-1) - radab(ji,jk) 
     640               WRITE(numout,*) ' dfc_i  : ', fc_i(ji,jk) - fc_i(ji,jk-1) 
     641               WRITE(numout,*) ' tot f  : ', fc_i(ji,jk) - fc_i(ji,jk-1) - radab(ji,jk) 
    660642            END DO 
    661643 
     
    667649 
    668650 
    669    SUBROUTINE lim_thd_con_dh(kideb,kiut,jl) 
     651   SUBROUTINE lim_thd_con_dh( kideb, kiut, jl ) 
    670652      !!----------------------------------------------------------------------- 
    671653      !!                   ***  ROUTINE lim_thd_con_dh  ***  
    672654      !!                  
    673655      !! ** Purpose :   Test energy conservation after enthalpy redistr. 
    674       !! 
    675       !! history : 
    676       !!  9.9  ! 07-04 (M.Vancoppenolle) original code 
    677656      !!----------------------------------------------------------------------- 
    678657      INTEGER, INTENT(in) ::        & 
     
    750729      ! Write ice state in case of big errors 
    751730      !--------------------------------------- 
    752  
    753731      DO ji = kideb, kiut 
    754732         IF ( cons_error(ji,jl) .GT. max_cons_err  ) THEN 
    755733            zji                 = MOD( npb(ji) - 1, jpi ) + 1 
    756734            zjj                 = ( npb(ji) - 1 ) / jpi + 1 
    757  
     735            ! 
    758736            WRITE(numout,*) ' alerte 1 - category : ', jl 
    759737            WRITE(numout,*) ' Untolerated conservation error after limthd_ent ' 
     
    776754            WRITE(numout,*) ' qt_s_in    : ', qt_s_in(ji,jl) / rdt_ice 
    777755            WRITE(numout,*) ' qt_i_in    : ', qt_i_in(ji,jl) / rdt_ice 
    778             WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + & 
    779                qt_s_in(ji,jl) ) / rdt_ice 
     756            WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) / rdt_ice 
    780757            WRITE(numout,*) ' qt_s_fin   : ', qt_s_fin(ji,jl) / rdt_ice 
    781758            WRITE(numout,*) ' qt_i_fin   : ', qt_i_fin(ji,jl) / rdt_ice 
    782             WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + & 
    783                qt_s_fin(ji,jl) ) / rdt_ice 
     759            WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) / rdt_ice 
    784760            WRITE(numout,*) ' * ' 
    785761            WRITE(numout,*) ' Ice variables --- : ' 
     
    790766            WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji) 
    791767            WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 
    792  
    793768         ENDIF 
    794769         ! 
     
    805780      !! 
    806781      !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
    807       !! 
    808782      !!------------------------------------------------------------------- 
    809783      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop 
     
    832806 
    833807   SUBROUTINE lim_thd_init 
    834  
    835808      !!----------------------------------------------------------------------- 
    836809      !!                   ***  ROUTINE lim_thd_init ***  
     
    843816      !! 
    844817      !! ** input   :   Namelist namicether 
    845        !!------------------------------------------------------------------- 
     818      !!------------------------------------------------------------------- 
    846819      NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb,   & 
    847820         &                hicmin, hiclim, amax  ,                                & 
Note: See TracChangeset for help on using the changeset viewer.