Changeset 10415


Ignore:
Timestamp:
2018-12-19T12:28:25+01:00 (22 months ago)
Author:
clem
Message:

correct the advection for ice ponds, make a couple of cosmetic changes and fix the conservation issue for sea ice advection (just a bugged control print)

Location:
NEMO/trunk/src
Files:
7 edited

Legend:

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

    r10413 r10415  
    332332   !! * Old values of global variables 
    333333   !!---------------------------------------------------------------------- 
    334    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b    !: snow and ice volumes/thickness 
    335    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b         !: 
    336    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                         !: snow heat content 
    337    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                         !: ice temperatures 
    338    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b              !: ice velocity 
    339    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                        !: ice concentration (total) 
     334   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b, h_ip_b    !: snow and ice volumes/thickness 
     335   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b                 !: 
     336   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                                 !: snow heat content 
     337   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                                 !: ice temperatures 
     338   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b                      !: ice velocity 
     339   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                                !: ice concentration (total) 
    340340             
    341341   !!---------------------------------------------------------------------- 
     
    441441      ! * Old values of global variables 
    442442      ii = ii + 1 
    443       ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl)        ,   & 
    444          &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,   & 
     443      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl), h_ip_b(jpi,jpj,jpl),  & 
     444         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,               & 
    445445         &      oa_i_b(jpi,jpj,jpl)                                                   , STAT=ierr(ii) ) 
    446446 
  • NEMO/trunk/src/ICE/icectl.F90

    r10069 r10415  
    190190 
    191191      ! heat flux 
    192       zhfx  = glob_sum( ( qt_atm_oi - qt_oce_ai - diag_heat - diag_trp_ei - diag_trp_es   & 
    193       !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check) 
    194          &              ) * e1e2t ) * zconv 
     192      zhfx  = glob_sum( ( qt_atm_oi - qt_oce_ai - diag_heat ) * e1e2t ) * zconv 
    195193 
    196194      ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
  • NEMO/trunk/src/ICE/icedyn.F90

    r10413 r10415  
    7676      INTEGER  ::   ji, jj, jl        ! dummy loop indices 
    7777      REAL(wp) ::   zcoefu, zcoefv 
    78       REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max 
     78      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
    7979      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i 
    8080      !!-------------------------------------------------------------------- 
     
    101101         DO jj = 2, jpjm1 
    102102            DO ji = fs_2, fs_jpim1 
    103                zhi_max(ji,jj,jl) = MAX( epsi20, h_i_b(ji,jj,jl), h_i_b(ji+1,jj  ,jl), h_i_b(ji  ,jj+1,jl), & 
    104                   &                                              h_i_b(ji-1,jj  ,jl), h_i_b(ji  ,jj-1,jl), & 
    105                   &                                              h_i_b(ji+1,jj+1,jl), h_i_b(ji-1,jj-1,jl), & 
    106                   &                                              h_i_b(ji+1,jj-1,jl), h_i_b(ji-1,jj+1,jl) ) 
    107                zhs_max(ji,jj,jl) = MAX( epsi20, h_s_b(ji,jj,jl), h_s_b(ji+1,jj  ,jl), h_s_b(ji  ,jj+1,jl), & 
    108                   &                                              h_s_b(ji-1,jj  ,jl), h_s_b(ji  ,jj-1,jl), & 
    109                   &                                              h_s_b(ji+1,jj+1,jl), h_s_b(ji-1,jj-1,jl), & 
    110                   &                                              h_s_b(ji+1,jj-1,jl), h_s_b(ji-1,jj+1,jl) ) 
     103               zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), & 
     104                  &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), & 
     105                  &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), & 
     106                  &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) ) 
     107               zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj  ,jl), h_i_b (ji  ,jj+1,jl), & 
     108                  &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), & 
     109                  &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), & 
     110                  &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) ) 
     111               zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj  ,jl), h_s_b (ji  ,jj+1,jl), & 
     112                  &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), & 
     113                  &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), & 
     114                  &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) ) 
    111115            END DO 
    112116         END DO 
    113117      END DO 
    114       CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1. ) 
     118      CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 
    115119      ! 
    116120      ! 
     
    118122 
    119123      CASE ( np_dynALL )           !==  all dynamical processes  ==! 
    120          CALL ice_dyn_rhg   ( kt )                                       ! -- rheology   
    121          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness 
    122          CALL ice_dyn_rdgrft( kt )                                       ! -- ridging/rafting  
    123          CALL ice_cor       ( kt , 1 )                                   ! -- Corrections 
     124         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
     125         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
     126         CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting  
     127         CALL ice_cor       ( kt , 1 )                                             ! -- Corrections 
    124128 
    125129      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==! 
    126          CALL ice_dyn_rhg   ( kt )                                       ! -- rheology   
    127          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness 
    128          CALL Hpiling                                                    ! -- simple pile-up (replaces ridging/rafting) 
     130         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
     131         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
     132         CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting) 
    129133 
    130134      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D) 
     
    163167         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 
    164168         ! --- 
    165          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness 
     169         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
    166170 
    167171         ! diagnostics: divergence at T points  
     
    186190 
    187191 
    188    SUBROUTINE Hbig( phi_max, phs_max ) 
     192   SUBROUTINE Hbig( phi_max, phs_max, phip_max ) 
    189193      !!------------------------------------------------------------------- 
    190194      !!                  ***  ROUTINE Hbig  *** 
     
    203207      !! ** input   : Max thickness of the surrounding 9-points 
    204208      !!------------------------------------------------------------------- 
    205       REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max   ! max ice thick from surrounding 9-pts 
     209      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    206210      ! 
    207211      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    208       REAL(wp) ::   zhi, zhs, zvs_excess, zfra 
     212      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra 
    209213      !!------------------------------------------------------------------- 
    210214      ! 
     
    216220               IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    217221                  ! 
     222                  !                               ! -- check h_ip -- ! 
     223                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
     224                  IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN 
     225                     zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) ) 
     226                     IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN 
     227                        a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     228                     ENDIF 
     229                  ENDIF 
     230                  ! 
    218231                  !                               ! -- check h_i -- ! 
    219232                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
    220233                  zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    221 !!clem                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)   
    222 !!clem                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. & 
    223 !!clem                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN 
    224234                  IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
    225235                     a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
     
    233243                     ! 
    234244                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice 
    235                      hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
     245                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
    236246                     ! 
    237                      e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra 
    238                      v_s(ji,jj,jl)   = a_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     247                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 
     248                     v_s(ji,jj,jl)          = a_i(ji,jj,jl) * phs_max(ji,jj,jl) 
    239249                  ENDIF            
    240250                  ! 
     
    242252                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 
    243253                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 
    244                   !    this imposed mini can artificially make the snow thickness very high (if concentration decreases drastically) 
     254                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 
    245255                  zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
    246256                  IF( zvs_excess > 0._wp ) THEN 
    247257                     zfra = zvs_excess / MAX( v_s(ji,jj,jl), epsi20 ) 
    248258                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 
    249                      hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
     259                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
    250260                     ! 
    251                      e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra 
    252                      v_s(ji,jj,jl)   = v_s(ji,jj,jl) - zvs_excess 
     261                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 
     262                     v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess 
    253263                  ENDIF 
    254264                   
  • NEMO/trunk/src/ICE/icedyn_rhg_evp.F90

    r10413 r10415  
    147147      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components 
    148148      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence 
    149       REAL(wp), DIMENSION(jpi,jpj) ::   zssh_lead_m                     ! array used for the calculation of ice surface slope: 
     149      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope: 
    150150      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
    151       !                                                                 !    ice top surface if ice is embedded    
     151      !                                                                 !    ice bottom surface if ice is embedded    
    152152      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array 
    153153      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array 
     
    268268      ! 2) Wind / ocean stress, mass terms, coriolis terms 
    269269      !------------------------------------------------------------------------------! 
    270  
    271       !== embedded sea ice: compute representative ice top surface      ==! 
    272       !== non-embedded sea ice: use ocean surface for slope calculation ==! 
    273       zssh_lead_m(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 
     270      ! sea surface height 
     271      !    embedded sea ice: compute representative ice top surface 
     272      !    non-embedded sea ice: use ocean surface for slope calculation 
     273      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 
    274274 
    275275      DO jj = 2, jpjm1 
     
    309309 
    310310            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
    311             zspgU(ji,jj)    = - zmassU * grav * ( zssh_lead_m(ji+1,jj) - zssh_lead_m(ji,jj) ) * r1_e1u(ji,jj) 
    312             zspgV(ji,jj)    = - zmassV * grav * ( zssh_lead_m(ji,jj+1) - zssh_lead_m(ji,jj) ) * r1_e2v(ji,jj) 
     311            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj) 
     312            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 
    313313 
    314314            ! masks 
  • NEMO/trunk/src/ICE/icestp.F90

    r10069 r10415  
    194194                                        CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling) 
    195195         ! 
    196 !! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work 
    197 !!       moreover it should only be called at the update frequency (as in agrif_ice_update.F90) 
    198 !# if defined key_agrif 
    199 !         IF( .NOT. Agrif_Root() )     CALL Agrif_ChildGrid_To_ParentGrid() 
    200 !# endif 
    201196                                        CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes 
    202 !# if defined key_agrif 
    203 !         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid() 
    204 !# endif 
     197         ! 
    205198         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs  
    206199         ! 
     
    368361      e_i_b (:,:,:,:) = e_i (:,:,:,:)   ! ice thermal energy 
    369362      WHERE( a_i_b(:,:,:) >= epsi20 ) 
    370          h_i_b(:,:,:) = v_i_b (:,:,:) / a_i_b(:,:,:)   ! ice thickness 
    371          h_s_b(:,:,:) = v_s_b (:,:,:) / a_i_b(:,:,:)   ! snw thickness 
     363         h_i_b(:,:,:) = v_i_b(:,:,:) / a_i_b(:,:,:)   ! ice thickness 
     364         h_s_b(:,:,:) = v_s_b(:,:,:) / a_i_b(:,:,:)   ! snw thickness 
    372365      ELSEWHERE 
    373366         h_i_b(:,:,:) = 0._wp 
    374367         h_s_b(:,:,:) = 0._wp 
     368      END WHERE 
     369       
     370      WHERE( a_ip(:,:,:) >= epsi20 ) 
     371         h_ip_b(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)   ! ice pond thickness 
     372      ELSEWHERE 
     373         h_ip_b(:,:,:) = 0._wp 
    375374      END WHERE 
    376375      ! 
  • NEMO/trunk/src/ICE/icevar.F90

    r10413 r10415  
    954954   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 
    955955      !!--------------------------------------------------------------------- 
    956       !!                   ***  ROUTINE rhg_evp_rst  *** 
     956      !!                   ***  ROUTINE ice_var_sshdyn  *** 
    957957      !!                      
    958958      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded 
     
    963963      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*,  
    964964      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008 
    965       !! 
    966965      !!---------------------------------------------------------------------- 
    967966      ! 
  • NEMO/trunk/src/OCE/stpctl.F90

    r10364 r10415  
    5050      !!              - Print it each 50 time steps 
    5151      !!              - Stop the run IF problem encountered by setting indic=-3 
    52       !!                Problems checked: |ssh| maximum larger than 10 m 
     52      !!                Problems checked: |ssh| maximum larger than 20 m 
    5353      !!                                  |U|   maximum larger than 10 m/s  
    5454      !!                                  negative sea surface salinity 
     
    136136      ENDIF 
    137137      ! 
    138       IF (  zmax(1) >   15._wp .OR.   &                    ! too large sea surface height ( > 15 m ) 
     138      IF (  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m ) 
    139139         &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
    140140         &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity 
     
    159159         IF(lwp) THEN 
    160160            WRITE(numout,cform_err) 
    161             WRITE(numout,*) ' stp_ctl: |ssh| > 10 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
     161            WRITE(numout,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
    162162            WRITE(numout,*) ' ======= ' 
    163163            WRITE(numout,9100) kt,   zmax(1), iih , ijh 
Note: See TracChangeset for help on using the changeset viewer.