New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 14037 for NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icestp.F90 – NEMO

Ignore:
Timestamp:
2020-12-03T12:20:38+01:00 (3 years ago)
Author:
ayoung
Message:

Updated to trunk at 14020. Sette tests passed with change of results for configurations with non-linear ssh. Ticket #2506.

Location:
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG

    • Property svn:externals

      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@13292        sette 
         10^/utils/CI/sette_wave@13990         sette 

  • NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icestp.F90

    r13216 r14037  
    5555   USE icedyn         ! sea-ice: dynamics 
    5656   USE icethd         ! sea-ice: thermodynamics 
    57    USE icecor         ! sea-ice: corrections 
    5857   USE iceupdate      ! sea-ice: sea surface boundary condition update 
    5958   USE icedia         ! sea-ice: budget diagnostics 
    … …  
    8685   PUBLIC   ice_init   ! called by sbcmod.F90 
    8786 
     87   !! * Substitutions 
     88#  include "do_loop_substitute.h90" 
    8889   !!---------------------------------------------------------------------- 
    8990   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    … …  
    120121      !!---------------------------------------------------------------------- 
    121122      ! 
    122       IF( ln_timing )   CALL timing_start('ice_stp') 
     123      IF( ln_timing )   CALL timing_start('icestp') 
    123124      ! 
    124125      !                                      !-----------------------! 
    … …  
    160161         IF( ln_icedyn .AND. .NOT.lk_c1d )   & 
    161162            &                           CALL ice_dyn( kt, Kmm )       ! -- Ice dynamics 
     163         ! 
     164                                        CALL diag_trends( 1 )         ! record dyn trends 
    162165         ! 
    163166         !                          !==  lateral boundary conditions  ==! 
    … …  
    188191         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics       
    189192         ! 
    190                                         CALL ice_cor( kt , 2 )        ! -- Corrections 
    191          ! 
     193                                        CALL diag_trends( 2 )         ! record thermo trends 
    192194                                        CALL ice_var_glo2eqv          ! necessary calls (at least for coupling) 
    193195                                        CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling) 
    … …  
    197199         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs  
    198200         ! 
     201         IF( ln_icediachk )             CALL ice_drift_wri( kt )      ! -- Diagnostics outputs for conservation  
     202         ! 
    199203                                        CALL ice_wri( kt )            ! -- Ice outputs  
    200204         ! 
    201205         IF( lrst_ice )                 CALL ice_rst_write( kt )      ! -- Ice restart file  
    202206         ! 
    203          IF( ln_icectl )                CALL ice_ctl( kt )            ! -- alerts in case of model crash 
     207         IF( ln_icectl )                CALL ice_ctl( kt )            ! -- Control checks 
    204208         ! 
    205209      ENDIF   ! End sea-ice time step only 
    … …  
    208212      ! --- Ocean time step --- ! 
    209213      !-------------------------! 
    210       IF( ln_icedyn )                   CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) )   ! -- update surface ocean stresses 
     214      CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) )         ! -- update surface ocean stresses 
    211215!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    212216      ! 
    213       IF( ln_timing )   CALL timing_stop('ice_stp') 
     217      IF( ln_timing )   CALL timing_stop('icestp') 
    214218      ! 
    215219   END SUBROUTINE ice_stp 
    … …  
    224228      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
    225229      ! 
    226       INTEGER :: ji, jj, ierr 
     230      INTEGER ::   ierr 
    227231      !!---------------------------------------------------------------------- 
    228232      IF(lwp) WRITE(numout,*) 
    … …  
    252256      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays') 
    253257      ! 
    254       CALL ice_itd_init                ! ice thickness distribution initialization 
    255       ! 
    256       CALL ice_thd_init                ! set ice thermodynics parameters (clem: important to call it first for melt ponds) 
    257       ! 
    258       !                                ! Initial sea-ice state 
    259       IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst 
    260          CALL ice_istate_init 
    261          CALL ice_istate( nit000, Kbb, Kmm, Kaa ) 
    262       ELSE                                    ! start from a restart file 
    263          CALL ice_rst_read( Kbb, Kmm, Kaa ) 
    264       ENDIF 
    265       CALL ice_var_glo2eqv 
    266       CALL ice_var_agg(1) 
    267       ! 
    268       CALL ice_sbc_init                ! set ice-ocean and ice-atm. coupling parameters 
    269       ! 
    270       CALL ice_dyn_init                ! set ice dynamics parameters 
    271       ! 
    272       CALL ice_update_init             ! ice surface boundary condition 
    273       ! 
    274       CALL ice_alb_init                ! ice surface albedo 
    275       ! 
    276       CALL ice_dia_init                ! initialization for diags 
    277       ! 
    278       fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction 
    279       tn_ice(:,:,:) = t_su(:,:,:)      ! initialisation of surface temp for coupled simu 
    280       ! 
    281258      !                                ! set max concentration in both hemispheres 
    282259      WHERE( gphit(:,:) > 0._wp )   ;   rn_amax_2d(:,:) = rn_amax_n  ! NH 
    283260      ELSEWHERE                     ;   rn_amax_2d(:,:) = rn_amax_s  ! SH 
    284261      END WHERE 
    285  
    286       IF( ln_rstart )   CALL iom_close( numrir )  ! close input ice restart file 
     262      ! 
     263      CALL diag_set0                   ! set diag of mass, heat and salt fluxes to 0: needed for Agrif child grids 
     264      ! 
     265      CALL ice_itd_init                ! ice thickness distribution initialization 
     266      ! 
     267      CALL ice_thd_init                ! set ice thermodynics parameters (clem: important to call it first for melt ponds) 
     268      ! 
     269      CALL ice_sbc_init                ! set ice-ocean and ice-atm. coupling parameters 
     270      ! 
     271      CALL ice_istate_init             ! Initial sea-ice state 
     272      IF ( ln_rstart .OR. nn_iceini_file == 2 ) THEN 
     273         CALL ice_rst_read( Kbb, Kmm, Kaa )         ! start from a restart file 
     274      ELSE 
     275         CALL ice_istate( nit000, Kbb, Kmm, Kaa )   ! start from rest or read a file 
     276      ENDIF 
     277      CALL ice_var_glo2eqv 
     278      CALL ice_var_agg(1) 
     279      ! 
     280      CALL ice_dyn_init                ! set ice dynamics parameters 
     281      ! 
     282      CALL ice_update_init             ! ice surface boundary condition 
     283      ! 
     284      CALL ice_alb_init                ! ice surface albedo 
     285      ! 
     286      CALL ice_dia_init                ! initialization for diags 
     287      ! 
     288      CALL ice_drift_init              ! initialization for diags of conservation 
     289      ! 
     290      fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction 
     291      tn_ice(:,:,:) = t_su(:,:,:)      ! initialisation of surface temp for coupled simu 
     292      ! 
     293      IF( ln_rstart )  THEN 
     294         CALL iom_close( numrir )  ! close input ice restart file 
     295         IF(lrxios) CALL iom_context_finalize(      cr_icerst_cxt         ) 
     296      ENDIF 
    287297      ! 
    288298   END SUBROUTINE ice_init 
    … …  
    340350      ENDIF 
    341351      ! 
    342       IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY') 
    343       ! 
    344352      rDt_ice   = REAL(nn_fsbc) * rn_Dt          !--- sea-ice timestep and its inverse 
    345353      r1_Dt_ice = 1._wp / rDt_ice 
    … …  
    365373      v_i_b (:,:,:)   = v_i (:,:,:)     ! ice volume 
    366374      v_s_b (:,:,:)   = v_s (:,:,:)     ! snow volume 
     375      v_ip_b(:,:,:)   = v_ip(:,:,:)     ! pond volume 
     376      v_il_b(:,:,:)   = v_il(:,:,:)     ! pond lid volume 
    367377      sv_i_b(:,:,:)   = sv_i(:,:,:)     ! salt content 
    368       oa_i_b(:,:,:)   = oa_i(:,:,:)     ! areal age content 
    369378      e_s_b (:,:,:,:) = e_s (:,:,:,:)   ! snow thermal energy 
    370379      e_i_b (:,:,:,:) = e_i (:,:,:,:)   ! ice thermal energy 
    … …  
    376385         h_s_b(:,:,:) = 0._wp 
    377386      END WHERE 
    378        
    379       WHERE( a_ip(:,:,:) >= epsi20 ) 
    380          h_ip_b(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)   ! ice pond thickness 
    381       ELSEWHERE 
    382          h_ip_b(:,:,:) = 0._wp 
    383       END WHERE 
    384387      ! 
    385388      ! ice velocities & total concentration 
    … …  
    398401      !!               of the time step 
    399402      !!---------------------------------------------------------------------- 
    400       INTEGER  ::   ji, jj      ! dummy loop index 
    401       !!---------------------------------------------------------------------- 
    402       sfx    (:,:) = 0._wp   ; 
    403       sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp 
    404       sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
    405       sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
    406       sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
    407       sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp 
    408       ! 
    409       wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
    410       wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
    411       wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
    412       wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
    413       wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
    414       wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp   
    415       wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp 
    416       wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp 
    417       wfx_snw_sni(:,:) = 0._wp  
    418       wfx_pnd(:,:) = 0._wp 
    419  
    420       hfx_thd(:,:) = 0._wp   ; 
    421       hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
    422       hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
    423       hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
    424       hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
    425       hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
    426       hfx_err_rem(:,:) = 0._wp 
    427       hfx_err_dif(:,:) = 0._wp 
    428       wfx_err_sub(:,:) = 0._wp 
    429       ! 
    430       diag_heat(:,:) = 0._wp ;   diag_sice(:,:) = 0._wp 
    431       diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp 
    432  
    433       ! SIMIP diagnostics 
    434       qcn_ice_bot(:,:,:) = 0._wp ; qcn_ice_top(:,:,:) = 0._wp ! conductive fluxes 
    435       t_si       (:,:,:) = rt0   ! temp at the ice-snow interface 
    436  
    437       tau_icebfr (:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
    438       cnd_ice    (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 
    439       qcn_ice    (:,:,:) = 0._wp   ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T) 
    440       qtr_ice_bot(:,:,:) = 0._wp   ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 
    441       qsb_ice_bot(:,:)   = 0._wp   ! (needed if ln_icethd=F) 
    442       ! 
    443       ! for control checks (ln_icediachk) 
    444       diag_trp_vi(:,:) = 0._wp   ;   diag_trp_vs(:,:) = 0._wp 
    445       diag_trp_ei(:,:) = 0._wp   ;   diag_trp_es(:,:) = 0._wp 
    446       diag_trp_sv(:,:) = 0._wp 
     403      INTEGER  ::   ji, jj, jl      ! dummy loop index 
     404      !!---------------------------------------------------------------------- 
     405 
     406      DO_2D( 1, 1, 1, 1 ) 
     407         sfx    (ji,jj) = 0._wp   ; 
     408         sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp 
     409         sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp 
     410         sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp 
     411         sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp 
     412         sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp 
     413         ! 
     414         wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp 
     415         wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp 
     416         wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp 
     417         wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp 
     418         wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp 
     419         wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp   
     420         wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp 
     421         wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp 
     422         wfx_snw_sni(ji,jj) = 0._wp  
     423         wfx_pnd(ji,jj) = 0._wp 
     424 
     425         hfx_thd(ji,jj) = 0._wp   ; 
     426         hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp 
     427         hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp 
     428         hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp 
     429         hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp 
     430         hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp 
     431         hfx_err_dif(ji,jj) = 0._wp 
     432         wfx_err_sub(ji,jj) = 0._wp 
     433         ! 
     434         diag_heat(ji,jj) = 0._wp ;   diag_sice(ji,jj) = 0._wp 
     435         diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp 
     436         diag_aice(ji,jj) = 0._wp ;   diag_vpnd(ji,jj) = 0._wp 
     437 
     438         tau_icebfr (ji,jj) = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
     439         qsb_ice_bot(ji,jj) = 0._wp   ! (needed if ln_icethd=F) 
     440 
     441         fhld(ji,jj) = 0._wp   ! needed if ln_icethd=F 
     442 
     443         ! for control checks (ln_icediachk) 
     444         diag_trp_vi(ji,jj) = 0._wp   ;   diag_trp_vs(ji,jj) = 0._wp 
     445         diag_trp_ei(ji,jj) = 0._wp   ;   diag_trp_es(ji,jj) = 0._wp 
     446         diag_trp_sv(ji,jj) = 0._wp 
     447         ! 
     448         diag_adv_mass(ji,jj) = 0._wp 
     449         diag_adv_salt(ji,jj) = 0._wp 
     450         diag_adv_heat(ji,jj) = 0._wp 
     451      END_2D 
     452 
     453      DO jl = 1, jpl 
     454         DO_2D( 1, 1, 1, 1 ) 
     455            ! SIMIP diagnostics 
     456            t_si       (ji,jj,jl) = rt0     ! temp at the ice-snow interface 
     457            qcn_ice_bot(ji,jj,jl) = 0._wp 
     458            qcn_ice_top(ji,jj,jl) = 0._wp   ! conductive fluxes 
     459            cnd_ice    (ji,jj,jl) = 0._wp   ! effective conductivity at the top of ice/snow (ln_cndflx=T) 
     460            qcn_ice    (ji,jj,jl) = 0._wp   ! conductive flux (ln_cndflx=T & ln_cndemule=T) 
     461            qtr_ice_bot(ji,jj,jl) = 0._wp   ! part of solar radiation transmitted through the ice needed at least for outputs 
     462            ! Melt pond surface melt diagnostics (mv - more efficient: grouped into one water volume flux) 
     463            dh_i_sum_2d(ji,jj,jl) = 0._wp 
     464            dh_s_mlt_2d(ji,jj,jl) = 0._wp 
     465         END_2D 
     466      ENDDO 
    447467 
    448468   END SUBROUTINE diag_set0 
     469 
     470 
     471   SUBROUTINE diag_trends( kn ) 
     472      !!---------------------------------------------------------------------- 
     473      !!                  ***  ROUTINE diag_trends  *** 
     474      !! 
     475      !! ** purpose : diagnostics of the trends. Used for conservation purposes 
     476      !!              and outputs 
     477      !!---------------------------------------------------------------------- 
     478      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo 
     479      !!---------------------------------------------------------------------- 
     480      ! 
     481      ! --- trends of heat, salt, mass (used for conservation controls) 
     482      IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN 
     483         ! 
     484         diag_heat(:,:) = diag_heat(:,:) & 
     485            &             - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & 
     486            &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 
     487         diag_sice(:,:) = diag_sice(:,:) & 
     488            &             + SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
     489         diag_vice(:,:) = diag_vice(:,:) & 
     490            &             + SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
     491         diag_vsnw(:,:) = diag_vsnw(:,:) & 
     492            &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos 
     493         diag_vpnd(:,:) = diag_vpnd(:,:) & 
     494            &             + SUM(     v_ip + v_il          - v_ip_b - v_il_b                , dim=3 ) * r1_Dt_ice * rhow 
     495         ! 
     496         IF( kn == 2 )    CALL iom_put ( 'hfxdhc' , diag_heat )   ! output of heat trend 
     497         ! 
     498      ENDIF 
     499      ! 
     500      ! --- trends of concentration (used for simip outputs) 
     501      IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 
     502         ! 
     503         diag_aice(:,:) = diag_aice(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 
     504         ! 
     505         IF( kn == 1 )   CALL iom_put( 'afxdyn' , diag_aice )                                           ! dyn trend 
     506         IF( kn == 2 )   CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice ) ! thermo trend 
     507         IF( kn == 2 )   CALL iom_put( 'afxtot' , diag_aice )                                           ! total trend 
     508         ! 
     509      ENDIF 
     510      ! 
     511   END SUBROUTINE diag_trends 
    449512 
    450513#else 
Note: See TracChangeset for help on using the changeset viewer.