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 13899 for NEMO/branches/2020/tickets_icb_1900/src/ICE/icestp.F90 – NEMO

Ignore:
Timestamp:
2020-11-27T17:26:33+01:00 (4 years ago)
Author:
mathiot
Message:

ticket #1900: update branch to trunk and add ICB test case

Location:
NEMO/branches/2020/tickets_icb_1900
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/tickets_icb_1900

    • Property svn:externals

      •  

        old new  
        22^/utils/build/makenemo@HEAD   makenemo 
        33^/utils/build/mk@HEAD         mk 
        4 ^/utils/tools/@HEAD           tools 
         4^/utils/tools@HEAD            tools 
        55^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        … …  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@12931        sette 
         10^/utils/CI/sette@13559        sette 

  • NEMO/branches/2020/tickets_icb_1900/src/ICE/icestp.F90

    r13216 r13899  
    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) 
    … …  
    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      ! 
    … …  
    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  
     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      ! 
    286293      IF( ln_rstart )   CALL iom_close( numrir )  ! close input ice restart file 
    287294      ! 
    … …  
    340347      ENDIF 
    341348      ! 
    342       IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY') 
    343       ! 
    344349      rDt_ice   = REAL(nn_fsbc) * rn_Dt          !--- sea-ice timestep and its inverse 
    345350      r1_Dt_ice = 1._wp / rDt_ice 
    … …  
    366371      v_s_b (:,:,:)   = v_s (:,:,:)     ! snow volume 
    367372      sv_i_b(:,:,:)   = sv_i(:,:,:)     ! salt content 
    368       oa_i_b(:,:,:)   = oa_i(:,:,:)     ! areal age content 
    369373      e_s_b (:,:,:,:) = e_s (:,:,:,:)   ! snow thermal energy 
    370374      e_i_b (:,:,:,:) = e_i (:,:,:,:)   ! ice thermal energy 
    … …  
    376380         h_s_b(:,:,:) = 0._wp 
    377381      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 
    384382      ! 
    385383      ! ice velocities & total concentration 
    … …  
    398396      !!               of the time step 
    399397      !!---------------------------------------------------------------------- 
    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 
     398      INTEGER  ::   ji, jj, jl      ! dummy loop index 
     399      !!---------------------------------------------------------------------- 
     400 
     401      DO_2D( 1, 1, 1, 1 ) 
     402         sfx    (ji,jj) = 0._wp   ; 
     403         sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp 
     404         sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp 
     405         sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp 
     406         sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp 
     407         sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp 
     408         ! 
     409         wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp 
     410         wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp 
     411         wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp 
     412         wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp 
     413         wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp 
     414         wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp   
     415         wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp 
     416         wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp 
     417         wfx_snw_sni(ji,jj) = 0._wp  
     418         wfx_pnd(ji,jj) = 0._wp 
     419 
     420         hfx_thd(ji,jj) = 0._wp   ; 
     421         hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp 
     422         hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp 
     423         hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp 
     424         hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp 
     425         hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp 
     426         hfx_err_dif(ji,jj) = 0._wp 
     427         wfx_err_sub(ji,jj) = 0._wp 
     428         ! 
     429         diag_heat(ji,jj) = 0._wp ;   diag_sice(ji,jj) = 0._wp 
     430         diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp 
     431 
     432         tau_icebfr (ji,jj) = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
     433         qsb_ice_bot(ji,jj) = 0._wp   ! (needed if ln_icethd=F) 
     434 
     435         fhld(ji,jj) = 0._wp   ! needed if ln_icethd=F 
     436 
     437         ! for control checks (ln_icediachk) 
     438         diag_trp_vi(ji,jj) = 0._wp   ;   diag_trp_vs(ji,jj) = 0._wp 
     439         diag_trp_ei(ji,jj) = 0._wp   ;   diag_trp_es(ji,jj) = 0._wp 
     440         diag_trp_sv(ji,jj) = 0._wp 
     441         ! 
     442         diag_adv_mass(ji,jj) = 0._wp 
     443         diag_adv_salt(ji,jj) = 0._wp 
     444         diag_adv_heat(ji,jj) = 0._wp 
     445      END_2D 
     446 
     447      DO jl = 1, jpl 
     448         DO_2D( 1, 1, 1, 1 ) 
     449            ! SIMIP diagnostics 
     450            t_si       (ji,jj,jl) = rt0     ! temp at the ice-snow interface 
     451            qcn_ice_bot(ji,jj,jl) = 0._wp 
     452            qcn_ice_top(ji,jj,jl) = 0._wp   ! conductive fluxes 
     453            cnd_ice    (ji,jj,jl) = 0._wp   ! effective conductivity at the top of ice/snow (ln_cndflx=T) 
     454            qcn_ice    (ji,jj,jl) = 0._wp   ! conductive flux (ln_cndflx=T & ln_cndemule=T) 
     455            qtr_ice_bot(ji,jj,jl) = 0._wp   ! part of solar radiation transmitted through the ice needed at least for outputs 
     456         END_2D 
     457      ENDDO 
    447458 
    448459   END SUBROUTINE diag_set0 
     460 
     461 
     462   SUBROUTINE diag_trends( kn ) 
     463      !!---------------------------------------------------------------------- 
     464      !!                  ***  ROUTINE diag_trends  *** 
     465      !! 
     466      !! ** purpose : diagnostics of the trends. Used for conservation purposes 
     467      !!              and outputs 
     468      !!---------------------------------------------------------------------- 
     469      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo 
     470      !!---------------------------------------------------------------------- 
     471      ! 
     472      ! --- trends of heat, salt, mass (used for conservation controls) 
     473      IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN 
     474         ! 
     475         diag_heat(:,:) = diag_heat(:,:) & 
     476            &             - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & 
     477            &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 
     478         diag_sice(:,:) = diag_sice(:,:) & 
     479            &             + SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
     480         diag_vice(:,:) = diag_vice(:,:) & 
     481            &             + SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
     482         diag_vsnw(:,:) = diag_vsnw(:,:) & 
     483            &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos 
     484         ! 
     485         IF( kn == 2 )    CALL iom_put ( 'hfxdhc' , diag_heat )   ! output of heat trend 
     486         ! 
     487      ENDIF 
     488      ! 
     489      ! --- trends of concentration (used for simip outputs) 
     490      IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 
     491         ! 
     492         diag_aice(:,:) = diag_aice(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 
     493         ! 
     494         IF( kn == 1 )   CALL iom_put( 'afxdyn' , diag_aice )                                           ! dyn trend 
     495         IF( kn == 2 )   CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice ) ! thermo trend 
     496         IF( kn == 2 )   CALL iom_put( 'afxtot' , diag_aice )                                           ! total trend 
     497         ! 
     498      ENDIF 
     499      ! 
     500   END SUBROUTINE diag_trends 
    449501 
    450502#else 
Note: See TracChangeset for help on using the changeset viewer.