Ignore:
Timestamp:
2017-09-22T16:55:24+02:00 (3 years ago)
Author:
clem
Message:

almost useless commits

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

    r8534 r8559  
    4242   !!   ice_var_eqv2glo   : transform from VEQV to VGLO 
    4343   !!   ice_var_salprof   : salinity profile in the ice 
    44    !!   ice_var_bv        : brine volume 
    4544   !!   ice_var_salprof1d : salinity profile in the ice 1D 
    4645   !!   ice_var_zapsmall  : remove very small area and volume 
    4746   !!   ice_var_itd       : convert 1-cat to multiple cat 
     47   !!   ice_var_bv        : brine volume 
    4848   !!---------------------------------------------------------------------- 
    4949   USE dom_oce        ! ocean space and time domain 
     
    6464   PUBLIC   ice_var_eqv2glo       
    6565   PUBLIC   ice_var_salprof       
    66    PUBLIC   ice_var_bv            
    6766   PUBLIC   ice_var_salprof1d     
    6867   PUBLIC   ice_var_zapsmall 
    6968   PUBLIC   ice_var_itd 
     69   PUBLIC   ice_var_bv            
    7070 
    7171   !!---------------------------------------------------------------------- 
     
    170170      END WHERE 
    171171      ! 
     172      WHERE( v_i(:,:,:) > epsi20 )   ;   z1_v_i(:,:,:) = 1._wp / v_i(:,:,:) 
     173      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp 
     174      END WHERE 
     175      ! 
    172176      ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness 
    173177 
     
    177181         ht_i  (:,:,jpl) = zhmax 
    178182         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax  
    179          z1_a_i(:,:,jpl) = zhmax / v_i(:,:,jpl)          ! NB: v_i always /=0 as ht_i > hi_max 
     183         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)          ! NB: v_i always /=0 as ht_i > hi_max 
    180184      END WHERE 
    181185 
     
    185189 
    186190      IF( nn_icesal == 2 ) THEN                    !--- salinity (with a minimum value imposed everywhere) 
    187          WHERE( v_i(:,:,:) > epsi20 )   ;   sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) ) 
     191         WHERE( v_i(:,:,:) > epsi20 )   ;   sm_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, smv_i(:,:,:) * z1_v_i(:,:,:) ) ) 
    188192         ELSEWHERE                      ;   sm_i(:,:,:) = rn_simin 
    189193         END WHERE 
     
    202206                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area  
    203207                     ! 
    204                      ze_i             =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3] 
     208                     ze_i             =   e_i(ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3] 
    205209                     ztmelts          = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C] 
    206210                     ! Conversion q(S,T) -> T (second order equation) 
     
    295299         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) ) 
    296300         ! 
    297          DO jk = 1, nlay_i 
    298             s_i(:,:,jk,:)  = sm_i(:,:,:) 
     301         DO jl = 1, jpl 
     302            DO jk = 1, nlay_i 
     303               s_i(:,:,jk,jl)  = sm_i(:,:,jl) 
     304            END DO 
    299305         END DO 
    300306         !                                      ! Slope of the linear profile  
     
    353359   END SUBROUTINE ice_var_salprof 
    354360 
    355  
    356    SUBROUTINE ice_var_bv 
    357       !!------------------------------------------------------------------- 
    358       !!                ***  ROUTINE ice_var_bv *** 
    359       !! 
    360       !! ** Purpose :   computes mean brine volume (%) in sea ice 
    361       !! 
    362       !! ** Method  : e = - 0.054 * S (ppt) / T (C) 
    363       !! 
    364       !! References : Vancoppenolle et al., JGR, 2007 
    365       !!------------------------------------------------------------------- 
    366       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    367       !!------------------------------------------------------------------- 
    368       ! 
    369 !!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done 
    370 !!   instead of setting everything to zero as just below 
    371       bv_i (:,:,:) = 0._wp 
    372       DO jl = 1, jpl 
    373          DO jk = 1, nlay_i 
    374             WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )    
    375                bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 
    376             END WHERE 
    377          END DO 
    378       END DO 
    379       WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 
    380       ELSEWHERE                     ;   bvm_i(:,:) = 0._wp 
    381       END WHERE 
    382       ! 
    383    END SUBROUTINE ice_var_bv 
    384  
    385  
    386361   SUBROUTINE ice_var_salprof1d 
    387362      !!------------------------------------------------------------------- 
     
    393368      INTEGER  ::   ji, jk    ! dummy loop indices 
    394369      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars 
    395       REAL(wp) ::   zalpha, zs, zs0              !   -      - 
    396       ! 
    397       REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s   ! 
     370      REAL(wp) ::   zs, zs0              !   -      - 
     371      ! 
     372      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   ! 
    398373      REAL(wp), PARAMETER :: zsi0 = 3.5_wp 
    399374      REAL(wp), PARAMETER :: zsi1 = 4.5_wp 
     
    405380      CASE( 1 )       !  constant salinity in time and space  ! 
    406381         !            !---------------------------------------! 
    407          s_i_1d(:,:) = rn_icesal 
     382         s_i_1d(1:nidx,:) = rn_icesal 
    408383         ! 
    409384         !            !---------------------------------------------! 
     
    411386         !            !---------------------------------------------! 
    412387         ! 
    413          ALLOCATE( z_slope_s(jpij) ) 
    414          ! 
    415          DO ji = 1, nidx          ! Slope of the linear profile 
    416             rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 
    417             z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) ) 
    418          END DO 
    419  
     388         ALLOCATE( z_slope_s(jpij), zalpha(jpij) ) 
     389         ! 
     390         !                                      ! Slope of the linear profile  
     391         WHERE( ht_i_1d(1:nidx) > epsi20 )   ;   z_slope_s(1:nidx) = 2._wp * sm_i_1d(1:nidx) / ht_i_1d(1:nidx) 
     392         ELSEWHERE                           ;   z_slope_s(1:nidx) = 0._wp 
     393         END WHERE 
     394          
    420395         z1_dS = 1._wp / ( zsi1 - zsi0 ) 
     396         DO ji = 1, nidx 
     397            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp  )  ) 
     398            !                             ! force a constant profile when SSS too low (Baltic Sea) 
     399            IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp 
     400         END DO 
     401         ! 
     402         ! Computation of the profile 
    421403         DO jk = 1, nlay_i 
    422404            DO ji = 1, nidx 
    423                zalpha = MAX(  0._wp , MIN(  ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp  )  ) 
    424                IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) )   zalpha = 0._wp               ! cst profile when SSS too low (Baltic Sea) 
    425                ! 
    426                zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i   ! linear profile with 0 surface value 
    427                zs  = zalpha * zs0 + ( 1._wp - zalpha ) * sm_i_1d(ji)                      ! weighting the profile 
    428                s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )                     ! bounding salinity 
    429             END DO 
    430          END DO 
    431          ! 
    432          DEALLOCATE( z_slope_s ) 
     405               !                          ! linear profile with 0 surface value 
     406               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 
     407               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * sm_i_1d(ji) 
     408               s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) ) 
     409            END DO 
     410         END DO 
     411         ! 
     412         DEALLOCATE( z_slope_s, zalpha ) 
    433413 
    434414         !            !-------------------------------------------! 
     
    436416         !            !-------------------------------------------!                                   (mean = 2.30) 
    437417         ! 
    438          sm_i_1d(:) = 2.30_wp 
     418         sm_i_1d(1:nidx) = 2.30_wp 
    439419         ! 
    440420!!gm cf remark in ice_var_salprof routine, CASE( 3 ) 
     
    459439      !!------------------------------------------------------------------- 
    460440      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
    461       REAL(wp) ::   zsal, zvi, zvs, zei, zes, zvp 
     441      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch 
    462442      !!------------------------------------------------------------------- 
    463443      ! 
     
    467447         ! Zap ice energy and use ocean heat to melt ice 
    468448         !----------------------------------------------------------------- 
     449         WHERE( a_i(:,:,jl) > epsi20 )   ;   ht_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl) 
     450         ELSEWHERE                       ;   ht_i(:,:,jl) = 0._wp 
     451         END WHERE 
     452         ! 
     453         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. ht_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp 
     454         ELSEWHERE                                                                            ;   zswitch(:,:) = 1._wp 
     455         END WHERE 
     456          
    469457         DO jk = 1, nlay_i 
    470458            DO jj = 1 , jpj 
    471459               DO ji = 1 , jpi 
    472 !!gm I think we can do better/faster  :  to be discussed 
    473                   rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
    474                   rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 
    475                   rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  & 
    476                      &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 
    477                   zei              = e_i(ji,jj,jk,jl) 
    478                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch 
    479                   t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch ) 
    480460                  ! update exchanges with ocean 
    481                   hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 
     461                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 
     462                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj) 
     463                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 
    482464               END DO 
    483465            END DO 
     
    486468         DO jj = 1 , jpj 
    487469            DO ji = 1 , jpi 
    488                rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
    489                rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 
    490                rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  & 
    491                   &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 
    492                zsal = smv_i(ji,jj,  jl) 
    493                zvi  = v_i  (ji,jj,  jl) 
    494                zvs  = v_s  (ji,jj,  jl) 
    495                zes  = e_s  (ji,jj,1,jl) 
    496                IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw )   zvp  = v_ip (ji,jj  ,jl) 
     470               ! update exchanges with ocean 
     471               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * smv_i(ji,jj,jl)   * rhoic * r1_rdtice 
     472               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i  (ji,jj,jl)   * rhoic * r1_rdtice 
     473               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s  (ji,jj,jl)   * rhosn * r1_rdtice 
     474               hfx_res(ji,jj)  = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s  (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0 
    497475               !----------------------------------------------------------------- 
    498476               ! Zap snow energy  
    499477               !----------------------------------------------------------------- 
    500                t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch ) 
    501                e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch 
     478               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 
     479               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj) 
    502480 
    503481               !----------------------------------------------------------------- 
    504482               ! zap ice and snow volume, add water and salt to ocean 
    505483               !----------------------------------------------------------------- 
    506                ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj) 
    507                a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch 
    508                v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch 
    509                v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch 
    510                t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch ) 
    511                oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch 
    512                smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch 
    513  
    514                ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch 
    515                ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch 
    516  
    517                ! MV MP 2016 
    518                IF ( ln_pnd ) THEN  
    519                   a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch 
    520                   v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch 
    521                   IF ( ln_pnd_fw )   wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_ip(ji,jj,jl)  - zvp  ) * rhofw * r1_rdtice 
    522                ENDIF 
    523                ! END MV MP 2016 
    524  
    525                ! update exchanges with ocean 
    526                sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
    527                wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
    528                wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
    529                hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 
    530             END DO 
    531          END DO 
     484               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj) 
     485               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * zswitch(ji,jj) 
     486               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * zswitch(ji,jj) 
     487               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * zswitch(ji,jj) 
     488               t_su (ji,jj,jl) = t_su (ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) ) 
     489               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * zswitch(ji,jj) 
     490               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * zswitch(ji,jj) 
     491 
     492               ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * zswitch(ji,jj) 
     493               ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * zswitch(ji,jj) 
     494 
     495            END DO 
     496         END DO 
     497 
     498         IF( ln_pnd ) THEN  
     499            DO jj = 1 , jpj 
     500               DO ji = 1 , jpi 
     501                  IF( ln_pnd_fw )   & 
     502                     &   wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice 
     503                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
     504                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
     505               END DO 
     506            END DO 
     507         ENDIF 
     508          
    532509      END DO  
    533510 
     
    688665    END SUBROUTINE ice_var_itd 
    689666 
     667 
     668    SUBROUTINE ice_var_bv 
     669      !!------------------------------------------------------------------- 
     670      !!                ***  ROUTINE ice_var_bv *** 
     671      !! 
     672      !! ** Purpose :   computes mean brine volume (%) in sea ice 
     673      !! 
     674      !! ** Method  : e = - 0.054 * S (ppt) / T (C) 
     675      !! 
     676      !! References : Vancoppenolle et al., JGR, 2007 
     677      !!------------------------------------------------------------------- 
     678      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
     679      !!------------------------------------------------------------------- 
     680      ! 
     681!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done 
     682!!   instead of setting everything to zero as just below 
     683      bv_i (:,:,:) = 0._wp 
     684      DO jl = 1, jpl 
     685         DO jk = 1, nlay_i 
     686            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )    
     687               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 
     688            END WHERE 
     689         END DO 
     690      END DO 
     691      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 
     692      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp 
     693      END WHERE 
     694      ! 
     695   END SUBROUTINE ice_var_bv 
     696 
     697 
    690698#else 
    691699   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.