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/limitd_me.F90

    r4147 r4161  
    55   !!====================================================================== 
    66   !! History :  LIM  ! 2006-02  (M. Vancoppenolle) Original code  
    7    !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & fsalt_rpo 
     7   !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_mec 
    88   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation 
    99   !!---------------------------------------------------------------------- 
     
    1212   !!   'key_lim3'                                      LIM-3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    USE par_oce        ! ocean parameters 
    15    USE dom_oce        ! ocean domain 
    16    USE phycst         ! physical constants (ocean directory)  
    17    USE sbc_oce        ! surface boundary condition: ocean fields 
    18    USE thd_ice        ! LIM thermodynamics 
    19    USE ice            ! LIM variables 
    20    USE par_ice        ! LIM parameters 
    21    USE dom_ice        ! LIM domain 
    22    USE limthd_lac     ! LIM 
    23    USE limvar         ! LIM 
    24    USE limcons        ! LIM 
    25    USE in_out_manager ! I/O manager 
    26    USE lbclnk         ! lateral boundary condition - MPP exchanges 
    27    USE lib_mpp        ! MPP library 
    28    USE wrk_nemo       ! work arrays 
    29    USE prtctl         ! Print control 
    30    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    31    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
     14   USE par_oce          ! ocean parameters 
     15   USE dom_oce          ! ocean domain 
     16   USE phycst           ! physical constants (ocean directory)  
     17   USE sbc_oce          ! surface boundary condition: ocean fields 
     18   USE thd_ice          ! LIM thermodynamics 
     19   USE ice              ! LIM variables 
     20   USE par_ice          ! LIM parameters 
     21   USE dom_ice          ! LIM domain 
     22   USE limthd_lac       ! LIM 
     23   USE limvar           ! LIM 
     24   USE limcons          ! LIM 
     25   USE in_out_manager   ! I/O manager 
     26   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     27   USE lib_mpp          ! MPP library 
     28   USE wrk_nemo         ! work arrays 
     29   USE prtctl           ! Print control 
     30  ! Check budget (Rousset) 
     31   USE iom              ! I/O manager 
     32   USE lib_fortran     ! glob_sum 
     33   USE limdiahsb 
     34   USE timing          ! Timing 
    3235 
    3336   IMPLICIT NONE 
     
    6265   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier 
    6366   REAL(wp), PARAMETER ::   kraft   = 2.0_wp    ! rafting multipliyer 
     67   REAL(wp), PARAMETER ::   kamax   = 1.0 
    6468 
    6569   REAL(wp) ::   Cp                             !  
     
    141145      REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2) 
    142146      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
     147      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     148      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
     149      ! mass and salt flux (clem) 
     150      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume... 
    143151      !!----------------------------------------------------------------------------- 
     152      IF( nn_timing == 1 )  CALL timing_start('limitd_me') 
    144153 
    145154      CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
     155 
     156      CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem 
    146157 
    147158      IF( numit == nstart  )   CALL lim_itd_me_init   ! Initialization (first time-step only) 
     
    151162         CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ') 
    152163      ENDIF 
     164 
     165      IF( ln_limdyn ) THEN          !   Start ridging and rafting   ! 
     166      ! ------------------------------- 
     167      !- check conservation (C Rousset) 
     168      IF (ln_limdiahsb) THEN 
     169         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     170         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     171         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 
     172         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 
     173      ENDIF 
     174      !- check conservation (C Rousset) 
     175      ! ------------------------------- 
     176 
     177      ! mass and salt flux init (clem) 
     178      zviold(:,:,:) = v_i(:,:,:) 
     179      zvsold(:,:,:) = v_s(:,:,:) 
     180      zsmvold(:,:,:) = smv_i(:,:,:) 
    153181 
    154182      !-----------------------------------------------------------------------------! 
     
    204232            ! to give asum = 1.0 after ridging. 
    205233 
    206             divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep 
     234            divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep 
    207235 
    208236            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) 
     
    288316         DO jj = 1, jpj 
    289317            DO ji = 1, jpi 
    290                IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN 
     318               IF (ABS(asum(ji,jj) - kamax ) .LT. epsi11) THEN 
    291319                  closing_net(ji,jj) = 0._wp 
    292320                  opning     (ji,jj) = 0._wp 
    293321               ELSE 
    294322                  iterate_ridging    = 1 
    295                   divu_adv   (ji,jj) = (1._wp - asum(ji,jj)) * r1_rdtice 
     323                  divu_adv   (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice 
    296324                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 
    297325                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) ) 
     
    330358         DO ji = 1, jpi 
    331359 
    332             IF( ABS( asum(ji,jj) - 1.0 ) > epsi11 )  asum_error = .true. 
     360            IF(ABS(asum(ji,jj) - kamax) > epsi11 ) asum_error = .true. 
    333361 
    334362            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice 
     
    349377      DO jj = 1, jpj 
    350378         DO ji = 1, jpi 
    351             IF( ABS( asum(ji,jj) - 1._wp )  >  epsi11 ) THEN   ! there is a bug 
     379            IF( ABS( asum(ji,jj) - kamax)  >  epsi11 ) THEN   ! there is a bug 
    352380               WRITE(numout,*) ' ' 
    353381               WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
     
    377405      CALL lim_var_glo2eqv 
    378406      CALL lim_itd_me_zapsmall 
     407 
     408      !-------------------------------- 
     409      ! Update mass/salt fluxes (clem) 
     410      !-------------------------------- 
     411      DO jl = 1, jpl 
     412         DO jj = 1, jpj  
     413            DO ji = 1, jpi 
     414               diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 
     415               rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic  
     416               rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn  
     417               sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice  
     418            END DO 
     419         END DO 
     420      END DO 
    379421 
    380422      !----------------- 
     
    425467      ENDIF 
    426468 
     469      ! ------------------------------- 
     470      !- check conservation (C Rousset) 
     471      IF (ln_limdiahsb) THEN 
     472         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
     473         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 
     474  
     475         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 
     476         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
     477 
     478         zchk_vmin = glob_min(v_i) 
     479         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     480         zchk_amin = glob_min(a_i) 
     481        
     482         IF(lwp) THEN 
     483            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * rday) 
     484            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday) 
     485            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3) 
     486            IF ( zchk_amax >  kamax+epsi10  ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax 
     487            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin 
     488         ENDIF 
     489      ENDIF 
     490      !- check conservation (C Rousset) 
     491      ! ------------------------------- 
     492 
    427493      !-------------------------! 
    428494      ! Back to initial values 
     
    448514 
    449515      ! heat content has to be corrected before ice volume 
    450       DO jl = 1, jpl 
    451          DO jk = 1, nlay_i 
    452             DO jj = 1, jpj 
    453                DO ji = 1, jpi 
    454                   IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 
    455                      ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 
    456                      old_e_i(ji,jj,jk,jl)   = d_e_i_trp(ji,jj,jk,jl) 
    457                      d_e_i_trp(ji,jj,jk,jl) = 0._wp 
    458                   ENDIF 
    459                END DO 
    460             END DO 
    461          END DO 
    462       END DO 
    463  
    464       DO jl = 1, jpl 
    465          DO jj = 1, jpj 
    466             DO ji = 1, jpi 
    467                IF(  old_v_i  (ji,jj,jl) < epsi06  .AND. & 
    468                     d_v_i_trp(ji,jj,jl) > epsi06    ) THEN 
    469                   old_v_i   (ji,jj,jl)   = d_v_i_trp(ji,jj,jl) 
    470                   d_v_i_trp (ji,jj,jl)   = 0._wp 
    471                   old_a_i   (ji,jj,jl)   = d_a_i_trp(ji,jj,jl) 
    472                   d_a_i_trp (ji,jj,jl)   = 0._wp 
    473                   old_v_s   (ji,jj,jl)   = d_v_s_trp(ji,jj,jl) 
    474                   d_v_s_trp (ji,jj,jl)   = 0._wp 
    475                   old_e_s   (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 
    476                   d_e_s_trp (ji,jj,1,jl) = 0._wp 
    477                   old_oa_i  (ji,jj,jl)   = d_oa_i_trp(ji,jj,jl) 
    478                   d_oa_i_trp(ji,jj,jl)   = 0._wp 
    479                   IF(  num_sal == 2  )   old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 
    480                   d_smv_i_trp(ji,jj,jl)  = 0._wp 
    481                ENDIF 
    482             END DO 
    483          END DO 
    484       END DO 
     516!clem@order 
     517!      DO jl = 1, jpl 
     518!         DO jk = 1, nlay_i 
     519!            DO jj = 1, jpj 
     520!               DO ji = 1, jpi 
     521!                  IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 
     522!                     ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 
     523!                     old_e_i(ji,jj,jk,jl)   = d_e_i_trp(ji,jj,jk,jl) 
     524!                     d_e_i_trp(ji,jj,jk,jl) = 0._wp 
     525!                  ENDIF 
     526!               END DO 
     527!            END DO 
     528!         END DO 
     529!      END DO 
     530! 
     531!      DO jl = 1, jpl 
     532!         DO jj = 1, jpj 
     533!            DO ji = 1, jpi 
     534!               IF(  old_v_i  (ji,jj,jl) < epsi06  .AND. & 
     535!                    d_v_i_trp(ji,jj,jl) > epsi06    ) THEN 
     536!                  old_v_i   (ji,jj,jl)   = d_v_i_trp(ji,jj,jl) 
     537!                  d_v_i_trp (ji,jj,jl)   = 0._wp 
     538!                  old_a_i   (ji,jj,jl)   = d_a_i_trp(ji,jj,jl) 
     539!                  d_a_i_trp (ji,jj,jl)   = 0._wp 
     540!                  old_v_s   (ji,jj,jl)   = d_v_s_trp(ji,jj,jl) 
     541!                  d_v_s_trp (ji,jj,jl)   = 0._wp 
     542!                  old_e_s   (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 
     543!                  d_e_s_trp (ji,jj,1,jl) = 0._wp 
     544!                  old_oa_i  (ji,jj,jl)   = d_oa_i_trp(ji,jj,jl) 
     545!                  d_oa_i_trp(ji,jj,jl)   = 0._wp 
     546!                  IF(  num_sal == 2  )   old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 
     547!                  d_smv_i_trp(ji,jj,jl)  = 0._wp 
     548!               ENDIF 
     549!            END DO 
     550!         END DO 
     551!      END DO 
     552!clem@order 
     553      ENDIF  ! ln_limdyn=.true. 
    485554      ! 
    486555      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    487556      ! 
     557      CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem 
     558      ! 
     559      IF( nn_timing == 1 )  CALL timing_stop('limitd_me') 
    488560   END SUBROUTINE lim_itd_me 
    489561 
     
    10861158            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 
    10871159 
    1088             IF (afrac(ji,jj) > 1.0 + epsi11) THEN  !riging 
     1160            IF (afrac(ji,jj) > kamax + epsi11) THEN  !riging 
    10891161               large_afrac = .true. 
    1090             ELSEIF (afrac(ji,jj) > 1.0) THEN  ! roundoff error 
    1091                afrac(ji,jj) = 1.0 
     1162            ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error 
     1163               afrac(ji,jj) = kamax 
    10921164            ENDIF 
    1093             IF (afrft(ji,jj) > 1.0 + epsi11) THEN !rafting 
     1165            IF (afrft(ji,jj) > kamax + epsi11) THEN !rafting 
    10941166               large_afrft = .true. 
    1095             ELSEIF (afrft(ji,jj) > 1.0) THEN  ! roundoff error 
    1096                afrft(ji,jj) = 1.0 
     1167            ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error 
     1168               afrft(ji,jj) = kamax 
    10971169            ENDIF 
    10981170 
     
    11371209             
    11381210            !                                                             ! excess of salt is flushed into the ocean 
    1139             sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 
    1140  
    1141             rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic / rau0   ! increase in ice volume du to seawater frozen in voids 
    1142              
     1211            !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 
     1212 
     1213            !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic    ! gurvan: increase in ice volume du to seawater frozen in voids              
     1214 
    11431215            !------------------------------------             
    11441216            ! 3.6 Increment ridging diagnostics 
     
    11501222            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 
    11511223            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 
    1152             diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 
     1224            !clem diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 
    11531225            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 
    11541226 
     
    12171289 
    12181290               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    1219                ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / nlay_i 
     1291               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 
    12201292 
    12211293               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
     
    12401312               ji = indxi(ij) 
    12411313               jj = indxj(ij) 
    1242                IF( afrac(ji,jj) > 1.0 + epsi11 ) THEN  
     1314               IF( afrac(ji,jj) > kamax + epsi11 ) THEN  
    12431315                  WRITE(numout,*) '' 
    12441316                  WRITE(numout,*) ' ardg > a_i' 
     
    12521324               ji = indxi(ij) 
    12531325               jj = indxj(ij) 
    1254                IF( afrft(ji,jj) > 1.0 + epsi11 ) THEN  
     1326               IF( afrft(ji,jj) > kamax + epsi11 ) THEN  
    12551327                  WRITE(numout,*) '' 
    12561328                  WRITE(numout,*) ' arft > a_i' 
Note: See TracChangeset for help on using the changeset viewer.