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 8518 – NEMO

Changeset 8518


Ignore:
Timestamp:
2017-09-13T18:46:56+02:00 (7 years ago)
Author:
clem
Message:

changes in style - part6 - commits of the day

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO
Files:
19 edited

Legend:

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

    r8517 r8518  
    170170    
    171171   !                                     !!** ice-dynamics namelist (namice_dyn) ** 
    172    LOGICAL,  PUBLIC ::   ln_dynFULL       !:  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction) 
    173    LOGICAL,  PUBLIC ::   ln_dynRHGADV     !:  dyn.: no ridge/raft & no corrections  (rheology + advection) 
    174    LOGICAL,  PUBLIC ::   ln_dynADV        !:  dyn.: only advection w prescribed vel.(rn_uvice + advection) 
    175172   REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice 
    176173   REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress 
     
    180177   REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction (landfast ice)  
    181178   ! 
    182    !                                     !!** ice-rdige/raft namelist (namice_rdgrft) ** 
    183    LOGICAL , PUBLIC ::   ln_str_H79       !: ice strength parameterization (Hibler79) 
    184    REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength, Hibler JPO79 
    185    REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength 
    186    LOGICAL , PUBLIC ::   ln_str_R75       !: ice strength parameterization (Rothrock75) 
    187    REAL(wp), PUBLIC ::   rn_perdg         !: ridging work divided by pot. energy change in ridging 
    188    ! 
    189179   !                                     !!** ice-rheology namelist (namice_rhg) ** 
    190    LOGICAL , PUBLIC ::   ln_rhg_EVP       !: EVP rheology 
    191180   REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9 
    192181   REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve 
    193182   INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling 
    194183   REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp)  
    195    ! 
    196    !                                     !!** ice-advection namelist (namice_adv) ** 
    197    LOGICAL , PUBLIC ::   ln_adv_Pra       !: Prather        advection scheme 
    198    LOGICAL , PUBLIC ::   ln_adv_UMx       !: Ultimate-Macho advection scheme 
    199    INTEGER , PUBLIC ::   nn_UMx           !: order of the UMx advection scheme    
    200184   ! 
    201185   !                                     !!** ice-thermodynamics namelist (namice_thd) ** 
     
    230214 
    231215   !                                     !!** ice-salinity namelist (namice_sal) ** 
    232    LOGICAL , PUBLIC ::   ln_icedS         !: activate gravity drainage and flushing (T) or not (F) 
    233216   INTEGER , PUBLIC ::   nn_icesal        !: salinity configuration used in the model 
    234217   !                                      ! 1 - constant salinity in both space and time 
     
    236219   !                                      ! 3 - salinity profile, constant in time 
    237220   REAL(wp), PUBLIC ::   rn_icesal        !: bulk salinity (ppt) in case of constant salinity 
    238    REAL(wp), PUBLIC ::   rn_sal_gd        !: restoring salinity for gravity drainage [PSU] 
    239    REAL(wp), PUBLIC ::   rn_time_gd       !: restoring time constant for gravity drainage (= 20 days) [s] 
    240    REAL(wp), PUBLIC ::   rn_sal_fl        !: restoring salinity for flushing [PSU] 
    241    REAL(wp), PUBLIC ::   rn_time_fl       !: restoring time constant for gravity drainage (= 10 days) [s] 
    242221   REAL(wp), PUBLIC ::   rn_simax         !: maximum ice salinity [PSU] 
    243222   REAL(wp), PUBLIC ::   rn_simin         !: minimum ice salinity [PSU] 
     
    329308   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice   [W.m-2] 
    330309   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt                          [W.m-2] 
    331    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion             [W.m-2] 
    332310   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2] 
    333311   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping             [W.m-2] 
     
    494472         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  & 
    495473         &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,  & 
    496          &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) ,     &  
     474         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) ,     &  
    497475         &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld   (jpi,jpj) ,                        & 
    498476         &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) ,     & 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90

    r8486 r8518  
    4747   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_opw_1d 
    4848   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_snw_1d 
    49    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_1d 
    5049   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_rem_1d 
    5150   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_dif_1d 
     
    178177         &      rn_amax_1d(jpij) ,                                         & 
    179178         &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) ,                      & 
    180          &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) ,   & 
     179         &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) ,                      & 
    181180         &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , hfx_out_1d(jpij), STAT=ierr(ii) ) 
    182181      ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90

    r8517 r8518  
    4141   INTEGER, PARAMETER ::   np_advPRA = 1   ! Prather scheme 
    4242   INTEGER, PARAMETER ::   np_advUMx = 2   ! Ultimate-Macho scheme 
    43  
     43   ! 
     44   ! ** namelist (namice_adv) ** 
     45   LOGICAL ::   ln_adv_Pra   ! Prather        advection scheme 
     46   LOGICAL ::   ln_adv_UMx   ! Ultimate-Macho advection scheme 
     47   INTEGER ::   nn_UMx       ! order of the UMx advection scheme    
     48   ! 
    4449   !! * Substitution 
    4550#  include "vectopt_loop_substitute.h90" 
     
    8590      CASE( np_advUMx )                ! ULTIMATE-MACHO scheme ! 
    8691         !                             !-----------------------! 
    87          CALL ice_adv_umx( kt, u_ice, v_ice,  & 
     92         CALL ice_adv_umx( nn_UMx, kt, u_ice, v_ice,  & 
    8893            &              ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    8994      !                                !-----------------------! 
     
    107112      IF( iom_use('deitrp') )   CALL iom_put( "deitrp" , diag_trp_ei         )         ! advected ice enthalpy (W/m2) 
    108113      IF( iom_use('destrp') )   CALL iom_put( "destrp" , diag_trp_es         )         ! advected snw enthalpy (W/m2) 
     114 
     115      IF( lrst_ice ) THEN                       !* write Prather fields in the restart file 
     116         IF( ln_adv_Pra )   CALL adv_pra_rst( 'WRITE', kt ) 
     117      ENDIF 
    109118                
    110119      ! controls 
     
    158167      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' ) 
    159168      ! 
     169      IF( ln_adv_Pra )   CALL adv_pra_rst( 'READ' )  !* read or initialize all required files 
     170      ! 
    160171   END SUBROUTINE ice_adv_init 
    161172 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90

    r8512 r8518  
    2020   USE lbclnk         ! lateral boundary condition - MPP exchanges 
    2121   USE in_out_manager ! I/O manager 
     22   USE iom            ! I/O library 
    2223   USE prtctl         ! Print control 
    2324   USE lib_mpp        ! MPP library 
     
    2829 
    2930   PUBLIC   ice_adv_prather   ! called by iceadv 
     31   PUBLIC   adv_pra_rst       ! called by iceadv 
    3032 
    3133   !! * Substitutions 
     
    592594   END SUBROUTINE adv_y 
    593595 
     596   SUBROUTINE adv_pra_rst( cdrw, kt ) 
     597      !!--------------------------------------------------------------------- 
     598      !!                   ***  ROUTINE adv_pra_rst  *** 
     599      !!                      
     600      !! ** Purpose :   Read or write RHG file in restart file 
     601      !! 
     602      !! ** Method  :   use of IOM library 
     603      !!---------------------------------------------------------------------- 
     604      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     605      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step 
     606      ! 
     607      INTEGER ::   jk, jl   ! dummy loop indices 
     608      INTEGER ::   iter     ! local integer 
     609      INTEGER ::   id1      ! local integer 
     610      CHARACTER(len=25) ::   znam 
     611      CHARACTER(len=2)  ::   zchar, zchar1 
     612      REAL(wp), DIMENSION(jpi,jpj) :: z2d 
     613      !!---------------------------------------------------------------------- 
     614      ! 
     615      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize 
     616         !                                   ! --------------- 
     617         IF( ln_rstart ) THEN                   !* Read the restart file 
     618            ! 
     619            id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. ) 
     620            ! 
     621            IF( id1 > 0 ) THEN      ! fields exist 
     622               DO jl = 1, jpl  
     623                  WRITE(zchar,'(I2.2)') jl 
     624                  znam = 'sxice'//'_htc'//zchar 
     625                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     626                  sxice(:,:,jl) = z2d(:,:) 
     627                  znam = 'syice'//'_htc'//zchar 
     628                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     629                  syice(:,:,jl) = z2d(:,:) 
     630                  znam = 'sxxice'//'_htc'//zchar 
     631                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     632                  sxxice(:,:,jl) = z2d(:,:) 
     633                  znam = 'syyice'//'_htc'//zchar 
     634                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     635                  syyice(:,:,jl) = z2d(:,:) 
     636                  znam = 'sxyice'//'_htc'//zchar 
     637                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     638                  sxyice(:,:,jl) = z2d(:,:) 
     639                  znam = 'sxsn'//'_htc'//zchar 
     640                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     641                  sxsn(:,:,jl) = z2d(:,:) 
     642                  znam = 'sysn'//'_htc'//zchar 
     643                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     644                  sysn(:,:,jl) = z2d(:,:) 
     645                  znam = 'sxxsn'//'_htc'//zchar 
     646                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     647                  sxxsn(:,:,jl) = z2d(:,:) 
     648                  znam = 'syysn'//'_htc'//zchar 
     649                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     650                  syysn(:,:,jl) = z2d(:,:) 
     651                  znam = 'sxysn'//'_htc'//zchar 
     652                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     653                  sxysn(:,:,jl) = z2d(:,:) 
     654                  znam = 'sxa'//'_htc'//zchar 
     655                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     656                  sxa(:,:,jl) = z2d(:,:) 
     657                  znam = 'sya'//'_htc'//zchar 
     658                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     659                  sya(:,:,jl) = z2d(:,:) 
     660                  znam = 'sxxa'//'_htc'//zchar 
     661                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     662                  sxxa(:,:,jl) = z2d(:,:) 
     663                  znam = 'syya'//'_htc'//zchar 
     664                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     665                  syya(:,:,jl) = z2d(:,:) 
     666                  znam = 'sxya'//'_htc'//zchar 
     667                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     668                  sxya(:,:,jl) = z2d(:,:) 
     669                  znam = 'sxc0'//'_htc'//zchar 
     670                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     671                  sxc0(:,:,jl) = z2d(:,:) 
     672                  znam = 'syc0'//'_htc'//zchar 
     673                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     674                  syc0(:,:,jl) = z2d(:,:) 
     675                  znam = 'sxxc0'//'_htc'//zchar 
     676                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     677                  sxxc0(:,:,jl) = z2d(:,:) 
     678                  znam = 'syyc0'//'_htc'//zchar 
     679                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     680                  syyc0(:,:,jl) = z2d(:,:) 
     681                  znam = 'sxyc0'//'_htc'//zchar 
     682                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     683                  sxyc0(:,:,jl) = z2d(:,:) 
     684                  znam = 'sxsal'//'_htc'//zchar 
     685                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     686                  sxsal(:,:,jl) = z2d(:,:) 
     687                  znam = 'sysal'//'_htc'//zchar 
     688                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     689                  sysal(:,:,jl) = z2d(:,:) 
     690                  znam = 'sxxsal'//'_htc'//zchar 
     691                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     692                  sxxsal(:,:,jl) = z2d(:,:) 
     693                  znam = 'syysal'//'_htc'//zchar 
     694                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     695                  syysal(:,:,jl) = z2d(:,:) 
     696                  znam = 'sxysal'//'_htc'//zchar 
     697                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     698                  sxysal(:,:,jl) = z2d(:,:) 
     699                  znam = 'sxage'//'_htc'//zchar 
     700                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     701                  sxage(:,:,jl) = z2d(:,:) 
     702                  znam = 'syage'//'_htc'//zchar 
     703                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     704                  syage(:,:,jl) = z2d(:,:) 
     705                  znam = 'sxxage'//'_htc'//zchar 
     706                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     707                  sxxage(:,:,jl) = z2d(:,:) 
     708                  znam = 'syyage'//'_htc'//zchar 
     709                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     710                  syyage(:,:,jl) = z2d(:,:) 
     711                  znam = 'sxyage'//'_htc'//zchar 
     712                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     713                  sxyage(:,:,jl)= z2d(:,:) 
     714               END DO 
     715               ! MV MP 2016 
     716               IF ( nn_pnd_scheme > 0 ) THEN 
     717                  DO jl = 1, jpl  
     718                     WRITE(zchar,'(I2.2)') jl 
     719                     znam = 'sxap'//'_htc'//zchar 
     720                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     721                     sxap(:,:,jl) = z2d(:,:) 
     722                     znam = 'syap'//'_htc'//zchar 
     723                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     724                     syap(:,:,jl) = z2d(:,:) 
     725                     znam = 'sxxap'//'_htc'//zchar 
     726                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     727                     sxxap(:,:,jl) = z2d(:,:) 
     728                     znam = 'syyap'//'_htc'//zchar 
     729                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     730                     syyap(:,:,jl) = z2d(:,:) 
     731                     znam = 'sxyap'//'_htc'//zchar 
     732                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     733                     sxyap(:,:,jl) = z2d(:,:) 
     734 
     735                     znam = 'sxvp'//'_htc'//zchar 
     736                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     737                     sxvp(:,:,jl) = z2d(:,:) 
     738                     znam = 'syvp'//'_htc'//zchar 
     739                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     740                     syvp(:,:,jl) = z2d(:,:) 
     741                     znam = 'sxxvp'//'_htc'//zchar 
     742                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     743                     sxxvp(:,:,jl) = z2d(:,:) 
     744                     znam = 'syyvp'//'_htc'//zchar 
     745                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     746                     syyvp(:,:,jl) = z2d(:,:) 
     747                     znam = 'sxyvp'//'_htc'//zchar 
     748                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     749                     sxyvp(:,:,jl) = z2d(:,:) 
     750                  END DO 
     751               ENDIF 
     752               ! END MV MP 2016 
     753 
     754               CALL iom_get( numrir, jpdom_autoglo, 'sxopw ' ,  sxopw  ) 
     755               CALL iom_get( numrir, jpdom_autoglo, 'syopw ' ,  syopw  ) 
     756               CALL iom_get( numrir, jpdom_autoglo, 'sxxopw' ,  sxxopw ) 
     757               CALL iom_get( numrir, jpdom_autoglo, 'syyopw' ,  syyopw ) 
     758               CALL iom_get( numrir, jpdom_autoglo, 'sxyopw' ,  sxyopw ) 
     759 
     760               DO jl = 1, jpl  
     761                  WRITE(zchar,'(I2.2)') jl 
     762                  DO jk = 1, nlay_i  
     763                     WRITE(zchar1,'(I2.2)') jk 
     764                     znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 
     765                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     766                     sxe(:,:,jk,jl) = z2d(:,:) 
     767                     znam = 'sye'//'_il'//zchar1//'_htc'//zchar 
     768                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     769                     sye(:,:,jk,jl) = z2d(:,:) 
     770                     znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 
     771                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     772                     sxxe(:,:,jk,jl) = z2d(:,:) 
     773                     znam = 'syye'//'_il'//zchar1//'_htc'//zchar 
     774                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     775                     syye(:,:,jk,jl) = z2d(:,:) 
     776                     znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 
     777                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     778                     sxye(:,:,jk,jl) = z2d(:,:) 
     779                  END DO 
     780               END DO 
     781               ! 
     782            ELSE                                     ! start rheology from rest 
     783               IF(lwp) WRITE(numout,*) '   ==>>   previous run without Prather, set moments to 0' 
     784               sxopw (:,:) = 0._wp   ;   sxice (:,:,:) = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:) = 0._wp 
     785               syopw (:,:) = 0._wp   ;   syice (:,:,:) = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:) = 0._wp 
     786               sxxopw(:,:) = 0._wp   ;   sxxice(:,:,:) = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:) = 0._wp 
     787               syyopw(:,:) = 0._wp   ;   syyice(:,:,:) = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:) = 0._wp 
     788               sxyopw(:,:) = 0._wp   ;   sxyice(:,:,:) = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:) = 0._wp 
     789               ! 
     790               sxc0  (:,:,:) = 0._wp   ;   sxe  (:,:,:,:) = 0._wp   ;   sxsal  (:,:,:) = 0._wp   ;   sxage  (:,:,:) = 0._wp 
     791               syc0  (:,:,:) = 0._wp   ;   sye  (:,:,:,:) = 0._wp   ;   sysal  (:,:,:) = 0._wp   ;   syage  (:,:,:) = 0._wp 
     792               sxxc0 (:,:,:) = 0._wp   ;   sxxe (:,:,:,:) = 0._wp   ;   sxxsal (:,:,:) = 0._wp   ;   sxxage (:,:,:) = 0._wp 
     793               syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp 
     794               sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp 
     795               ! MV MP 2016 
     796               IF ( nn_pnd_scheme > 0 ) THEN 
     797                  sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp  
     798                  syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp  
     799                  sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp  
     800                  syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp  
     801                  sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp 
     802               ENDIF 
     803               ! END MV MP 2016 
     804            ENDIF 
     805         ELSE                                   !* Start from rest 
     806            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set moments to 0' 
     807            sxopw (:,:) = 0._wp   ;   sxice (:,:,:) = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:) = 0._wp 
     808            syopw (:,:) = 0._wp   ;   syice (:,:,:) = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:) = 0._wp 
     809            sxxopw(:,:) = 0._wp   ;   sxxice(:,:,:) = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:) = 0._wp 
     810            syyopw(:,:) = 0._wp   ;   syyice(:,:,:) = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:) = 0._wp 
     811            sxyopw(:,:) = 0._wp   ;   sxyice(:,:,:) = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:) = 0._wp 
     812            ! 
     813            sxc0  (:,:,:) = 0._wp   ;   sxe  (:,:,:,:) = 0._wp   ;   sxsal  (:,:,:) = 0._wp   ;   sxage  (:,:,:) = 0._wp 
     814            syc0  (:,:,:) = 0._wp   ;   sye  (:,:,:,:) = 0._wp   ;   sysal  (:,:,:) = 0._wp   ;   syage  (:,:,:) = 0._wp 
     815            sxxc0 (:,:,:) = 0._wp   ;   sxxe (:,:,:,:) = 0._wp   ;   sxxsal (:,:,:) = 0._wp   ;   sxxage (:,:,:) = 0._wp 
     816            syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp 
     817            sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp 
     818            ! MV MP 2016 
     819            IF ( nn_pnd_scheme > 0 ) THEN 
     820               sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp  
     821               syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp  
     822               sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp  
     823               syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp  
     824               sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp 
     825            ENDIF 
     826            ! END MV MP 2016 
     827         ENDIF 
     828         ! 
     829      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
     830         !                                   ! ------------------- 
     831         IF(lwp) WRITE(numout,*) '---- adv-rst ----' 
     832         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1 
     833         ! 
     834         DO jl = 1, jpl  
     835            WRITE(zchar,'(I2.2)') jl 
     836            znam = 'sxice'//'_htc'//zchar 
     837            z2d(:,:) = sxice(:,:,jl) 
     838            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     839            znam = 'syice'//'_htc'//zchar 
     840            z2d(:,:) = syice(:,:,jl) 
     841            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     842            znam = 'sxxice'//'_htc'//zchar 
     843            z2d(:,:) = sxxice(:,:,jl) 
     844            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     845            znam = 'syyice'//'_htc'//zchar 
     846            z2d(:,:) = syyice(:,:,jl) 
     847            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     848            znam = 'sxyice'//'_htc'//zchar 
     849            z2d(:,:) = sxyice(:,:,jl) 
     850            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     851            znam = 'sxsn'//'_htc'//zchar 
     852            z2d(:,:) = sxsn(:,:,jl) 
     853            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     854            znam = 'sysn'//'_htc'//zchar 
     855            z2d(:,:) = sysn(:,:,jl) 
     856            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     857            znam = 'sxxsn'//'_htc'//zchar 
     858            z2d(:,:) = sxxsn(:,:,jl) 
     859            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     860            znam = 'syysn'//'_htc'//zchar 
     861            z2d(:,:) = syysn(:,:,jl) 
     862            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     863            znam = 'sxysn'//'_htc'//zchar 
     864            z2d(:,:) = sxysn(:,:,jl) 
     865            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     866            znam = 'sxa'//'_htc'//zchar 
     867            z2d(:,:) = sxa(:,:,jl) 
     868            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     869            znam = 'sya'//'_htc'//zchar 
     870            z2d(:,:) = sya(:,:,jl) 
     871            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     872            znam = 'sxxa'//'_htc'//zchar 
     873            z2d(:,:) = sxxa(:,:,jl) 
     874            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     875            znam = 'syya'//'_htc'//zchar 
     876            z2d(:,:) = syya(:,:,jl) 
     877            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     878            znam = 'sxya'//'_htc'//zchar 
     879            z2d(:,:) = sxya(:,:,jl) 
     880            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     881            znam = 'sxc0'//'_htc'//zchar 
     882            z2d(:,:) = sxc0(:,:,jl) 
     883            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     884            znam = 'syc0'//'_htc'//zchar 
     885            z2d(:,:) = syc0(:,:,jl) 
     886            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     887            znam = 'sxxc0'//'_htc'//zchar 
     888            z2d(:,:) = sxxc0(:,:,jl) 
     889            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     890            znam = 'syyc0'//'_htc'//zchar 
     891            z2d(:,:) = syyc0(:,:,jl) 
     892            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     893            znam = 'sxyc0'//'_htc'//zchar 
     894            z2d(:,:) = sxyc0(:,:,jl) 
     895            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     896            znam = 'sxsal'//'_htc'//zchar 
     897            z2d(:,:) = sxsal(:,:,jl) 
     898            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     899            znam = 'sysal'//'_htc'//zchar 
     900            z2d(:,:) = sysal(:,:,jl) 
     901            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     902            znam = 'sxxsal'//'_htc'//zchar 
     903            z2d(:,:) = sxxsal(:,:,jl) 
     904            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     905            znam = 'syysal'//'_htc'//zchar 
     906            z2d(:,:) = syysal(:,:,jl) 
     907            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     908            znam = 'sxysal'//'_htc'//zchar 
     909            z2d(:,:) = sxysal(:,:,jl) 
     910            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     911            znam = 'sxage'//'_htc'//zchar 
     912            z2d(:,:) = sxage(:,:,jl) 
     913            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     914            znam = 'syage'//'_htc'//zchar 
     915            z2d(:,:) = syage(:,:,jl) 
     916            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     917            znam = 'sxxage'//'_htc'//zchar 
     918            z2d(:,:) = sxxage(:,:,jl) 
     919            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     920            znam = 'syyage'//'_htc'//zchar 
     921            z2d(:,:) = syyage(:,:,jl) 
     922            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     923            znam = 'sxyage'//'_htc'//zchar 
     924            z2d(:,:) = sxyage(:,:,jl) 
     925            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     926         END DO 
     927 
     928         CALL iom_rstput( iter, nitrst, numriw, 'sxopw ' ,  sxopw  ) 
     929         CALL iom_rstput( iter, nitrst, numriw, 'syopw ' ,  syopw  ) 
     930         CALL iom_rstput( iter, nitrst, numriw, 'sxxopw' ,  sxxopw ) 
     931         CALL iom_rstput( iter, nitrst, numriw, 'syyopw' ,  syyopw ) 
     932         CALL iom_rstput( iter, nitrst, numriw, 'sxyopw' ,  sxyopw ) 
     933          
     934         DO jl = 1, jpl  
     935            WRITE(zchar,'(I2.2)') jl 
     936            DO jk = 1, nlay_i  
     937               WRITE(zchar1,'(I2.2)') jk 
     938               znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 
     939               z2d(:,:) = sxe(:,:,jk,jl) 
     940               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     941               znam = 'sye'//'_il'//zchar1//'_htc'//zchar 
     942               z2d(:,:) = sye(:,:,jk,jl) 
     943               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     944               znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 
     945               z2d(:,:) = sxxe(:,:,jk,jl) 
     946               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     947               znam = 'syye'//'_il'//zchar1//'_htc'//zchar 
     948               z2d(:,:) = syye(:,:,jk,jl) 
     949               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     950               znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 
     951               z2d(:,:) = sxye(:,:,jk,jl) 
     952               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     953            END DO 
     954         END DO 
     955         ! MV MP 2016 
     956         IF ( nn_pnd_scheme > 0 ) THEN 
     957            DO jl = 1, jpl  
     958               WRITE(zchar,'(I2.2)') jl 
     959               znam = 'sxap'//'_htc'//zchar 
     960               z2d(:,:) = sxap(:,:,jl) 
     961               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     962               znam = 'syap'//'_htc'//zchar 
     963               z2d(:,:) = syap(:,:,jl) 
     964               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     965               znam = 'sxxap'//'_htc'//zchar 
     966               z2d(:,:) = sxxap(:,:,jl) 
     967               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     968               znam = 'syyap'//'_htc'//zchar 
     969               z2d(:,:) = syyap(:,:,jl) 
     970               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     971               znam = 'sxyap'//'_htc'//zchar 
     972               z2d(:,:) = sxyap(:,:,jl) 
     973               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     974    
     975               znam = 'sxvp'//'_htc'//zchar 
     976               z2d(:,:) = sxvp(:,:,jl) 
     977               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     978               znam = 'syvp'//'_htc'//zchar 
     979               z2d(:,:) = syvp(:,:,jl) 
     980               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     981               znam = 'sxxvp'//'_htc'//zchar 
     982               z2d(:,:) = sxxvp(:,:,jl) 
     983               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     984               znam = 'syyvp'//'_htc'//zchar 
     985               z2d(:,:) = syyvp(:,:,jl) 
     986               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     987               znam = 'sxyvp'//'_htc'//zchar 
     988               z2d(:,:) = sxyvp(:,:,jl) 
     989               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     990            END DO 
     991         ENDIF 
     992         ! 
     993      ENDIF 
     994      ! 
     995   END SUBROUTINE adv_pra_rst 
     996 
    594997#else 
    595998   !!---------------------------------------------------------------------- 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90

    r8512 r8518  
    2929   PRIVATE 
    3030 
    31    PUBLIC   ice_adv_umx    ! routine called by iceadv.F90 
     31   PUBLIC   ice_adv_umx    ! called by iceadv.F90 
    3232       
    3333   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6 
     
    4343CONTAINS 
    4444 
    45    SUBROUTINE ice_adv_umx( kt, pu_ice, pv_ice,  & 
     45   SUBROUTINE ice_adv_umx( k_order, kt, pu_ice, pv_ice,  & 
    4646      &                    pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    4747      !!---------------------------------------------------------------------- 
     
    5454      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    5555      !!---------------------------------------------------------------------- 
     56      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the scheme (1-5 or 20) 
    5657      INTEGER                     , INTENT(in   ) ::   kt         ! time step 
    5758      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
     
    112113      !---------------! 
    113114      DO jt = 1, initad 
    114          CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area  
     115         CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area  
    115116         DO jl = 1, jpl 
    116             CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area 
    117             CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume 
    118             CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) )       ! Salt content 
    119             CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) )       ! Age content 
     117            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area 
     118            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume 
     119            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) )       ! Salt content 
     120            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) )       ! Age content 
    120121            DO jk = 1, nlay_i 
    121                CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content 
    122             END DO 
    123             CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume 
    124             CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content 
     122               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content 
     123            END DO 
     124            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume 
     125            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content 
    125126            IF ( nn_pnd_scheme > 0 ) THEN 
    126                CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction 
    127                CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume 
     127               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction 
     128               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume 
    128129            ENDIF 
    129130         END DO 
     
    134135   END SUBROUTINE ice_adv_umx 
    135136    
    136    SUBROUTINE adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc ) 
     137   SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc ) 
    137138      !!---------------------------------------------------------------------- 
    138139      !!                  ***  ROUTINE adv_umx  *** 
     
    147148      !! ** Action : - pt  the after advective tracer 
    148149      !!---------------------------------------------------------------------- 
     150      INTEGER                     , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme 
    149151      INTEGER                     , INTENT(in   ) ::   kt             ! number of iteration 
    150152      REAL(wp)                    , INTENT(in   ) ::   pdt            ! tracer time-step 
     
    189191      ! High order (_ho) fluxes  
    190192      ! ----------------------- 
    191       SELECT CASE( nn_UMx ) 
     193      SELECT CASE( k_order ) 
    192194      CASE ( 20 )                          ! centered second order 
    193195         DO jj = 2, jpjm1 
     
    199201         ! 
    200202      CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme 
    201          CALL macho( kt, nn_UMx, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 
     203         CALL macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 
    202204         ! 
    203205         DO jj = 2, jpjm1 
     
    240242 
    241243 
    242    SUBROUTINE macho( kt, k_order, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 
     244   SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 
    243245      !!--------------------------------------------------------------------- 
    244246      !!                    ***  ROUTINE ultimate_x  *** 
     
    251253      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    252254      !!---------------------------------------------------------------------- 
     255      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme 
    253256      INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration 
    254       INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme 
    255257      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step 
    256258      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc        ! tracer fields 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90

    r8514 r8518  
    4545   !!---------------------------------------------------------------------- 
    4646CONTAINS 
     47 
     48   INTEGER FUNCTION ice_dia_alloc() 
     49      !!---------------------------------------------------------------------! 
     50      !!                ***  ROUTINE ice_rdgrft_alloc *** 
     51      !!---------------------------------------------------------------------! 
     52      ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ice_dia_alloc ) 
     53 
     54      IF( lk_mpp             )   CALL mpp_sum ( ice_dia_alloc ) 
     55      IF( ice_dia_alloc /= 0 )   CALL ctl_warn( 'ice_dia_alloc: failed to allocate arrays' ) 
     56      ! 
     57   END FUNCTION ice_dia_alloc 
    4758 
    4859   SUBROUTINE ice_dia( kt ) 
     
    188199      !       
    189200      IF( ln_icediahsb ) THEN 
    190          ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror ) 
    191          IF( ierror > 0 )  THEN 
    192             CALL ctl_stop( 'ice_dia: unable to allocate vol_loc_ini' ) 
    193             RETURN 
    194          ENDIF 
    195          ! 
    196          CALL ice_dia_rst( 'READ' )  !* read or initialize all required files 
     201         IF( ice_dia_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dia_init : unable to allocate arrays' )   ! allocate tke arrays 
     202         CALL ice_dia_rst( 'READ' )   ! read or initialize all required files 
    197203      ENDIF 
    198204      ! 
     
    220226            CALL iom_get( numrir, 'kt_ice' , ziter ) 
    221227            IF(lwp) WRITE(numout,*) 
    222             IF(lwp) WRITE(numout,*) ' ice_dia_rst read at time step = ', ziter 
    223             IF(lwp) WRITE(numout,*) '~~~~~~~' 
     228            IF(lwp) WRITE(numout,*) 'ice_dia_rst read at time step = ', ziter 
     229            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
    224230            CALL iom_get( numrir, 'frc_voltop' , frc_voltop  ) 
    225231            CALL iom_get( numrir, 'frc_volbot' , frc_volbot  ) 
     
    252258         IF( iter == nitrst ) THEN 
    253259            IF(lwp) WRITE(numout,*) 
    254             IF(lwp) WRITE(numout,*) ' ice_dia_rst write at time step = ', kt 
    255             IF(lwp) WRITE(numout,*) '~~~~~~~' 
     260            IF(lwp) WRITE(numout,*) 'ice_dia_rst write at time step = ', kt 
     261            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    256262         ENDIF 
    257263         ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90

    r8517 r8518  
    4242   ! 
    4343   ! ** namelist (namdyn) ** 
    44    REAL(wp) ::   rn_uice          ! prescribed u-vel (case np_dynADV) 
    45    REAL(wp) ::   rn_vice          ! prescribed v-vel (case np_dynADV) 
     44   LOGICAL  ::   ln_dynFULL       ! full ice dynamics               (rheology + advection + ridging/rafting + correction) 
     45   LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections  (rheology + advection) 
     46   LOGICAL  ::   ln_dynADV        ! only advection w prescribed vel.(rn_uvice + advection) 
     47   REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV) 
     48   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV) 
    4649    
    4750   !! * Substitutions 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90

    r8515 r8518  
    1717   !!   ice_istate_init  :  initialization of ice state and namelist read 
    1818   !!---------------------------------------------------------------------- 
     19   USE par_oce        ! ocean parameters 
    1920   USE phycst         ! physical constant 
    2021   USE oce            ! dynamics and tracers variables 
    2122   USE dom_oce        ! ocean domain 
    22    USE sbc_oce , ONLY : sst_m, sss_m  
    23    USE sbc_ice , ONLY : tn_ice 
     23   USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd  
     24   USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b 
    2425   USE eosbn2         ! equation of state 
     26   USE domvvl         ! Variable volume 
    2527   USE ice            ! sea-ice variables 
    26    USE par_oce        ! ocean parameters 
    2728   USE icevar         ! ice_var_salprof 
    2829   ! 
     
    9192      !!              where there is no ice (clem: I do not know why, is it mandatory?)  
    9293      !!-------------------------------------------------------------------- 
    93       INTEGER    :: ji, jj, jk, jl             ! dummy loop indices 
    94       REAL(wp)   :: ztmelts, zdh 
    95       INTEGER    :: i_hemis, i_fill, jl0 
    96       REAL(wp)   :: zarg, zV, zconv, zdv 
    97       REAL(wp), DIMENSION(jpi,jpj)     :: zswitch    ! ice indicator 
    98       REAL(wp), DIMENSION(jpi,jpj)     :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file 
    99       REAL(wp), DIMENSION(jpi,jpj)     :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 
    100       REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini, za_i_ini                         !data by cattegories to fill 
    101       INTEGER , DIMENSION(4)           :: itest 
     94      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices 
     95      REAL(wp) ::   ztmelts, zdh 
     96      INTEGER  ::   i_hemis, i_fill, jl0 
     97      REAL(wp) ::   zarg, zV, zconv, zdv 
     98      INTEGER , DIMENSION(4)           ::   itest 
     99      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d 
     100      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator 
     101      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file 
     102      REAL(wp), DIMENSION(jpi,jpj)     ::   zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 
     103      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini, za_i_ini                         !data by cattegories to fill 
    102104      !-------------------------------------------------------------------- 
    103105 
     
    119121      CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 
    120122      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)  
    121  
    122123 
    123124      !-------------------------------------------------------------------- 
     
    377378         ! Melt pond volume and fraction 
    378379         IF ( ln_pnd ) THEN 
    379  
    380380            DO jl = 1, jpl 
    381  
    382381               a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:) 
    383382               h_ip     (:,:,jl) = 0.05 * zswitch(:,:) 
    384383               a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl)  
    385384               v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl) 
    386  
    387             END DO 
    388  
     385            END DO 
    389386         ELSE 
    390  
    391387            a_ip(:,:,:)      = 0._wp 
    392388            v_ip(:,:,:)      = 0._wp 
    393389            a_ip_frac(:,:,:) = 0._wp 
    394390            h_ip     (:,:,:) = 0._wp 
    395  
    396391         ENDIF 
    397392         ! END MV MP 2016 
     
    437432      u_ice (:,:)     = 0._wp 
    438433      v_ice (:,:)     = 0._wp 
    439       stress1_i(:,:)  = 0._wp 
    440       stress2_i(:,:)  = 0._wp 
    441       stress12_i(:,:) = 0._wp 
    442  
    443       !-------------------------------------------------------------------- 
    444       ! 5) Moments for advection 
    445       !-------------------------------------------------------------------- 
    446  
    447       sxopw (:,:) = 0._wp  
    448       syopw (:,:) = 0._wp  
    449       sxxopw(:,:) = 0._wp  
    450       syyopw(:,:) = 0._wp  
    451       sxyopw(:,:) = 0._wp 
    452  
    453       sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp 
    454       syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp 
    455       sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp 
    456       syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp 
    457       sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp 
    458  
    459       sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp    
    460       syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp    
    461       sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp    
    462       syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp    
    463       sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp    
    464  
    465       sxsal  (:,:,:)  = 0._wp 
    466       sysal  (:,:,:)  = 0._wp 
    467       sxxsal (:,:,:)  = 0._wp 
    468       syysal (:,:,:)  = 0._wp 
    469       sxysal (:,:,:)  = 0._wp 
    470  
    471       sxage  (:,:,:)  = 0._wp 
    472       syage  (:,:,:)  = 0._wp 
    473       sxxage (:,:,:)  = 0._wp 
    474       syyage (:,:,:)  = 0._wp 
    475       sxyage (:,:,:)  = 0._wp 
    476  
     434      ! 
     435      !-------------------------------------------------------------------- 
     436      ! Snow-ice mass (case ice is fully embedded)                    |  
     437      !-------------------------------------------------------------------- 
     438      snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )   ! snow+ice mass 
     439      snwice_mass_b(:,:) = snwice_mass(:,:) 
     440      ! 
     441      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area 
     442 
     443         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
     444         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 
     445 
     446         IF( .NOT.ln_linssh ) THEN 
     447             
     448            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:) 
     449            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE 
     450          
     451            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors                 
     452               e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:) 
     453               e3t_b(:,:,jk) = e3t_n(:,:,jk) 
     454               e3t_a(:,:,jk) = e3t_n(:,:,jk) 
     455            END DO 
     456             
     457            ! Reconstruction of all vertical scale factors at now and before time-steps 
     458            ! ========================================================================= 
     459            ! Horizontal scale factor interpolations 
     460            ! -------------------------------------- 
     461            CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 
     462            CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 
     463            CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
     464            CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
     465            CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 
     466            ! Vertical scale factor interpolations 
     467            ! ------------------------------------ 
     468            CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  ) 
     469            CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
     470            CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
     471            CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
     472            CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
     473            ! t- and w- points depth 
     474            ! ---------------------- 
     475            !!gm not sure of that.... 
     476            gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
     477            gdepw_n(:,:,1) = 0.0_wp 
     478            gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
     479            DO jk = 2, jpk 
     480               gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  ) 
     481               gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) 
     482               gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:) 
     483            END DO 
     484         ENDIF 
     485      ENDIF 
     486       
    477487      !------------------------------------ 
    478488      ! 6) store fields at before time-step 
     
    488498      u_ice_b(:,:)     = u_ice(:,:) 
    489499      v_ice_b(:,:)     = v_ice(:,:) 
    490  
    491       ! MV MP 2016 
    492       IF ( nn_pnd_scheme > 0 ) THEN 
    493          sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp  
    494          syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp  
    495          sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp  
    496          syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp  
    497          sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp 
    498       ENDIF 
    499       ! END MV MP 2016 
    500500 
    501501!!!clem 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90

    r8517 r8518  
    3434   PUBLIC   ice_rdgrft_strength      ! called by icerhg_evp 
    3535   PUBLIC   ice_rdgrft_init          ! called by ice_stp 
    36    PUBLIC   ice_rdgrft_alloc         ! called by ice_init  
    3736 
    3837   ! Variables shared among ridging subroutines 
    39    REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
    40    REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
    41    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/closing associated w/ category n 
    42    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
    43    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
    44    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice 
    45    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! thickness of ridging ice / mean ridge thickness 
    46    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging 
    47    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting 
     38   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
     39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
     40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/closing associated w/ category n 
     41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
     42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
     43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice 
     44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   krdg     ! thickness of ridging ice / mean ridge thickness 
     45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting 
    4847   ! 
    4948   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier 
     
    5251   ! 
    5352   ! ** namelist (namice_rdgrft) ** 
     53   LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79) 
     54   REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79 
     55   REAL(wp) ::   rn_crhg          ! determines changes in ice strength 
     56   LOGICAL  ::   ln_str_R75       ! ice strength parameterization (Rothrock75) 
     57   REAL(wp) ::   rn_perdg         ! ridging work divided by pot. energy change in ridging 
    5458   REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging             
    5559   LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975)) 
     
    8286         &      hrmin(jpi,jpj,jpl) , hraft (jpi,jpj,jpl)   , aridge(jpi,jpj,jpl) ,     & 
    8387         &      hrmax(jpi,jpj,jpl) , krdg  (jpi,jpj,jpl)   , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc ) 
    84          ! 
     88 
     89      IF( lk_mpp                )   CALL mpp_sum ( ice_rdgrft_alloc ) 
    8590      IF( ice_rdgrft_alloc /= 0 )   CALL ctl_warn( 'ice_rdgrft_alloc: failed to allocate arrays' ) 
    8691      ! 
     
    909914         CALL ctl_stop( 'ice_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 
    910915      ENDIF 
     916      !                              ! allocate tke arrays 
     917      IF( ice_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_rdgrft_init : unable to allocate arrays' ) 
    911918      ! 
    912919  END SUBROUTINE ice_rdgrft_init 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90

    r8517 r8518  
    3737   INTEGER, PARAMETER ::   np_rhgEVP = 1   ! EVP rheology 
    3838!! INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology 
    39     
     39 
     40   ! ** namelist (namrhg) ** 
     41   LOGICAL ::   ln_rhg_EVP       ! EVP rheology 
     42   ! 
    4043   !! * Substitutions 
    4144#  include "vectopt_loop_substitute.h90" 
     
    7780         !                             !------------------------! 
    7881         CALL ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 
    79  
     82         !          
    8083      END SELECT 
     84      ! 
     85      IF( lrst_ice ) THEN                       !* write EVP fields in the restart file 
     86         IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt ) 
     87      ENDIF 
    8188      ! 
    8289      ! controls 
     
    8693      ! 
    8794   END SUBROUTINE ice_rhg 
    88  
    8995 
    9096   SUBROUTINE ice_rhg_init 
     
    131137!!    IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF 
    132138      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_rhg_init: choose one and only one ice rheology' ) 
    133  
    134      !                              ! allocate tke arrays 
    135 !!clem example      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' ) 
     139      ! 
     140      IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files 
    136141      ! 
    137142   END SUBROUTINE ice_rhg_init 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90

    r8517 r8518  
    4141   PRIVATE 
    4242 
    43    PUBLIC   ice_rhg_evp        ! routine called by icerhg.F90 
     43   PUBLIC   ice_rhg_evp   ! called by icerhg.F90 
     44   PUBLIC   rhg_evp_rst   ! called by icerhg.F90 
    4445 
    4546   !! * Substitutions 
     
    108109      !!------------------------------------------------------------------- 
    109110      INTEGER, INTENT(in) ::   kt      ! time step 
    110       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i 
    111       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i  
     111      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i 
     112      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i  
    112113      !! 
    113114      INTEGER ::   ji, jj       ! dummy loop indices 
     
    825826   END SUBROUTINE ice_rhg_evp 
    826827 
     828   SUBROUTINE rhg_evp_rst( cdrw, kt ) 
     829      !!--------------------------------------------------------------------- 
     830      !!                   ***  ROUTINE rhg_evp_rst  *** 
     831      !!                      
     832      !! ** Purpose :   Read or write RHG file in restart file 
     833      !! 
     834      !! ** Method  :   use of IOM library 
     835      !!---------------------------------------------------------------------- 
     836      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     837      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step 
     838      ! 
     839      INTEGER  ::   iter            ! local integer 
     840      INTEGER  ::   id1, id2, id3   ! local integers 
     841      !!---------------------------------------------------------------------- 
     842      ! 
     843      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize 
     844         !                                   ! --------------- 
     845         IF( ln_rstart ) THEN                   !* Read the restart file 
     846            ! 
     847            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. ) 
     848            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. ) 
     849            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. ) 
     850            ! 
     851            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist 
     852               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  ) 
     853               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  ) 
     854               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i ) 
     855            ELSE                                     ! start rheology from rest 
     856               IF(lwp) WRITE(numout,*) '   ==>>   previous run without rheology, set stresses to 0' 
     857               stress1_i (:,:) = 0._wp 
     858               stress2_i (:,:) = 0._wp 
     859               stress12_i(:,:) = 0._wp 
     860            ENDIF 
     861         ELSE                                   !* Start from rest 
     862            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set stresses to 0' 
     863            stress1_i (:,:) = 0._wp 
     864            stress2_i (:,:) = 0._wp 
     865            stress12_i(:,:) = 0._wp 
     866         ENDIF 
     867         ! 
     868      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
     869         !                                   ! ------------------- 
     870         IF(lwp) WRITE(numout,*) '---- rhg-rst ----' 
     871         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1 
     872         ! 
     873         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  ) 
     874         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  ) 
     875         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i ) 
     876         ! 
     877      ENDIF 
     878      ! 
     879   END SUBROUTINE rhg_evp_rst 
     880 
    827881#else 
    828882   !!---------------------------------------------------------------------- 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90

    r8514 r8518  
    1717   !!---------------------------------------------------------------------- 
    1818   USE ice            ! sea-ice variables 
    19    USE sbc_ice , ONLY :  snwice_mass, snwice_mass_b 
    2019   USE dom_oce        ! ocean domain 
    2120   USE sbc_oce , ONLY : nn_fsbc 
     
    3029   PRIVATE 
    3130 
    32    PUBLIC   ice_rst_opn    ! routine called by icestep.F90 
    33    PUBLIC   ice_rst_write  ! routine called by icestep.F90 
    34    PUBLIC   ice_rst_read   ! routine called by ice_init 
    35  
    36    LOGICAL, PUBLIC ::   lrst_ice         !: logical to control the ice restart write  
    37    INTEGER, PUBLIC ::   numrir, numriw   !: logical unit for ice restart (read and write) 
     31   PUBLIC   ice_rst_opn     ! called by icestp 
     32   PUBLIC   ice_rst_write   ! called by icestp 
     33   PUBLIC   ice_rst_read    ! called by ice_init 
    3834 
    3935   !!---------------------------------------------------------------------- 
     
    7975                  WRITE(numout,*) '             open ice restart NetCDF file: ',TRIM(clpath)//clname 
    8076               END SELECT 
    81                IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN    
    82                   WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp 
    83                ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp 
     77               IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN 
     78                  WRITE(numout,*) '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp 
     79               ELSE 
     80                  WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp 
    8481               ENDIF 
    8582            ENDIF 
     
    10299      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    103100      !! 
    104       INTEGER ::   ji, jj, jk ,jl   ! dummy loop indices 
     101      INTEGER ::   jk ,jl   ! dummy loop indices 
    105102      INTEGER ::   iter 
    106103      CHARACTER(len=25) ::   znam 
     
    127124!!gm   "just" ask Sebatien 
    128125 
    129  
    130126      ! Prognostic variables  
    131127      DO jl = 1, jpl  
     
    133129         znam = 'v_i'//'_htc'//zchar 
    134130         z2d(:,:) = v_i(:,:,jl) 
    135          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     131         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_i 
    136132         znam = 'v_s'//'_htc'//zchar 
    137133         z2d(:,:) = v_s(:,:,jl) 
    138          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     134         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_s 
    139135         znam = 'smv_i'//'_htc'//zchar 
    140136         z2d(:,:) = smv_i(:,:,jl) 
    141          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     137         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! smv_i 
    142138         znam = 'oa_i'//'_htc'//zchar 
    143139         z2d(:,:) = oa_i(:,:,jl) 
    144          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     140         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! oa_i 
    145141         znam = 'a_i'//'_htc'//zchar 
    146142         z2d(:,:) = a_i(:,:,jl) 
    147          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     143         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_i 
    148144         znam = 't_su'//'_htc'//zchar 
    149145         z2d(:,:) = t_su(:,:,jl) 
    150          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     146         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! t_su 
    151147      END DO 
    152148 
     
    157153            znam = 'a_ip'//'_htc'//zchar 
    158154            z2d(:,:) = a_ip(:,:,jl) 
    159             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     155            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_ip 
    160156            znam = 'v_ip'//'_htc'//zchar 
    161157            z2d(:,:) = v_ip(:,:,jl) 
    162             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     158            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_ip 
    163159         END DO 
    164160      ENDIF 
     
    169165         znam = 'tempt_sl1'//'_htc'//zchar 
    170166         z2d(:,:) = e_s(:,:,1,jl) 
    171          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     167         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! e_s 
    172168         DO jk = 1, nlay_i  
    173169            WRITE(zchar1,'(I2.2)') jk 
    174170            znam = 'tempt'//'_il'//zchar1//'_htc'//zchar 
    175171            z2d(:,:) = e_i(:,:,jk,jl) 
    176             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    177          END DO 
    178       END DO 
    179  
    180       CALL iom_rstput( iter, nitrst, numriw, 'u_ice'        , u_ice      ) 
    181       CALL iom_rstput( iter, nitrst, numriw, 'v_ice'        , v_ice      ) 
    182       CALL iom_rstput( iter, nitrst, numriw, 'stress1_i'    , stress1_i  ) 
    183       CALL iom_rstput( iter, nitrst, numriw, 'stress2_i'    , stress2_i  ) 
    184       CALL iom_rstput( iter, nitrst, numriw, 'stress12_i'   , stress12_i ) 
    185       CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass ) 
    186       CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) 
    187  
    188       ! In case Prather scheme is used for advection, write second order moments 
    189       ! ------------------------------------------------------------------------ 
    190       IF( ln_adv_Pra ) THEN 
    191           
    192          DO jl = 1, jpl  
    193             WRITE(zchar,'(I2.2)') jl 
    194             znam = 'sxice'//'_htc'//zchar 
    195             z2d(:,:) = sxice(:,:,jl) 
    196             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    197             znam = 'syice'//'_htc'//zchar 
    198             z2d(:,:) = syice(:,:,jl) 
    199             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    200             znam = 'sxxice'//'_htc'//zchar 
    201             z2d(:,:) = sxxice(:,:,jl) 
    202             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    203             znam = 'syyice'//'_htc'//zchar 
    204             z2d(:,:) = syyice(:,:,jl) 
    205             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    206             znam = 'sxyice'//'_htc'//zchar 
    207             z2d(:,:) = sxyice(:,:,jl) 
    208             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    209             znam = 'sxsn'//'_htc'//zchar 
    210             z2d(:,:) = sxsn(:,:,jl) 
    211             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    212             znam = 'sysn'//'_htc'//zchar 
    213             z2d(:,:) = sysn(:,:,jl) 
    214             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    215             znam = 'sxxsn'//'_htc'//zchar 
    216             z2d(:,:) = sxxsn(:,:,jl) 
    217             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    218             znam = 'syysn'//'_htc'//zchar 
    219             z2d(:,:) = syysn(:,:,jl) 
    220             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    221             znam = 'sxysn'//'_htc'//zchar 
    222             z2d(:,:) = sxysn(:,:,jl) 
    223             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    224             znam = 'sxa'//'_htc'//zchar 
    225             z2d(:,:) = sxa(:,:,jl) 
    226             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    227             znam = 'sya'//'_htc'//zchar 
    228             z2d(:,:) = sya(:,:,jl) 
    229             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    230             znam = 'sxxa'//'_htc'//zchar 
    231             z2d(:,:) = sxxa(:,:,jl) 
    232             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    233             znam = 'syya'//'_htc'//zchar 
    234             z2d(:,:) = syya(:,:,jl) 
    235             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    236             znam = 'sxya'//'_htc'//zchar 
    237             z2d(:,:) = sxya(:,:,jl) 
    238             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    239             znam = 'sxc0'//'_htc'//zchar 
    240             z2d(:,:) = sxc0(:,:,jl) 
    241             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    242             znam = 'syc0'//'_htc'//zchar 
    243             z2d(:,:) = syc0(:,:,jl) 
    244             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    245             znam = 'sxxc0'//'_htc'//zchar 
    246             z2d(:,:) = sxxc0(:,:,jl) 
    247             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    248             znam = 'syyc0'//'_htc'//zchar 
    249             z2d(:,:) = syyc0(:,:,jl) 
    250             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    251             znam = 'sxyc0'//'_htc'//zchar 
    252             z2d(:,:) = sxyc0(:,:,jl) 
    253             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    254             znam = 'sxsal'//'_htc'//zchar 
    255             z2d(:,:) = sxsal(:,:,jl) 
    256             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    257             znam = 'sysal'//'_htc'//zchar 
    258             z2d(:,:) = sysal(:,:,jl) 
    259             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    260             znam = 'sxxsal'//'_htc'//zchar 
    261             z2d(:,:) = sxxsal(:,:,jl) 
    262             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    263             znam = 'syysal'//'_htc'//zchar 
    264             z2d(:,:) = syysal(:,:,jl) 
    265             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    266             znam = 'sxysal'//'_htc'//zchar 
    267             z2d(:,:) = sxysal(:,:,jl) 
    268             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    269             znam = 'sxage'//'_htc'//zchar 
    270             z2d(:,:) = sxage(:,:,jl) 
    271             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    272             znam = 'syage'//'_htc'//zchar 
    273             z2d(:,:) = syage(:,:,jl) 
    274             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    275             znam = 'sxxage'//'_htc'//zchar 
    276             z2d(:,:) = sxxage(:,:,jl) 
    277             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    278             znam = 'syyage'//'_htc'//zchar 
    279             z2d(:,:) = syyage(:,:,jl) 
    280             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    281             znam = 'sxyage'//'_htc'//zchar 
    282             z2d(:,:) = sxyage(:,:,jl) 
    283             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    284          END DO 
    285  
    286          CALL iom_rstput( iter, nitrst, numriw, 'sxopw ' ,  sxopw  ) 
    287          CALL iom_rstput( iter, nitrst, numriw, 'syopw ' ,  syopw  ) 
    288          CALL iom_rstput( iter, nitrst, numriw, 'sxxopw' ,  sxxopw ) 
    289          CALL iom_rstput( iter, nitrst, numriw, 'syyopw' ,  syyopw ) 
    290          CALL iom_rstput( iter, nitrst, numriw, 'sxyopw' ,  sxyopw ) 
    291           
    292          DO jl = 1, jpl  
    293             WRITE(zchar,'(I2.2)') jl 
    294             DO jk = 1, nlay_i  
    295                WRITE(zchar1,'(I2.2)') jk 
    296                znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 
    297                z2d(:,:) = sxe(:,:,jk,jl) 
    298                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    299                znam = 'sye'//'_il'//zchar1//'_htc'//zchar 
    300                z2d(:,:) = sye(:,:,jk,jl) 
    301                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    302                znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 
    303                z2d(:,:) = sxxe(:,:,jk,jl) 
    304                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    305                znam = 'syye'//'_il'//zchar1//'_htc'//zchar 
    306                z2d(:,:) = syye(:,:,jk,jl) 
    307                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    308                znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 
    309                z2d(:,:) = sxye(:,:,jk,jl) 
    310                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    311             END DO 
    312          END DO 
    313          ! MV MP 2016 
    314          IF ( nn_pnd_scheme > 0 ) THEN 
    315             DO jl = 1, jpl  
    316                WRITE(zchar,'(I2.2)') jl 
    317                znam = 'sxap'//'_htc'//zchar 
    318                z2d(:,:) = sxap(:,:,jl) 
    319                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    320                znam = 'syap'//'_htc'//zchar 
    321                z2d(:,:) = syap(:,:,jl) 
    322                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    323                znam = 'sxxap'//'_htc'//zchar 
    324                z2d(:,:) = sxxap(:,:,jl) 
    325                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    326                znam = 'syyap'//'_htc'//zchar 
    327                z2d(:,:) = syyap(:,:,jl) 
    328                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    329                znam = 'sxyap'//'_htc'//zchar 
    330                z2d(:,:) = sxyap(:,:,jl) 
    331                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    332     
    333                znam = 'sxvp'//'_htc'//zchar 
    334                z2d(:,:) = sxvp(:,:,jl) 
    335                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    336                znam = 'syvp'//'_htc'//zchar 
    337                z2d(:,:) = syvp(:,:,jl) 
    338                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    339                znam = 'sxxvp'//'_htc'//zchar 
    340                z2d(:,:) = sxxvp(:,:,jl) 
    341                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    342                znam = 'syyvp'//'_htc'//zchar 
    343                z2d(:,:) = syyvp(:,:,jl) 
    344                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    345                znam = 'sxyvp'//'_htc'//zchar 
    346                z2d(:,:) = sxyvp(:,:,jl) 
    347                CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    348             END DO 
    349          ENDIF 
    350  
    351       ENDIF 
     172            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! e_i 
     173         END DO 
     174      END DO 
     175 
     176      CALL iom_rstput( iter, nitrst, numriw, 'u_ice', u_ice ) ! u_ice 
     177      CALL iom_rstput( iter, nitrst, numriw, 'v_ice', v_ice ) ! v_ice 
    352178       
    353179      ! close restart file 
     
    368194      !! ** purpose  :   read of sea-ice variable restart in a netcdf file 
    369195      !!---------------------------------------------------------------------- 
    370       INTEGER  :: ji, jj, jk, jl 
     196      INTEGER  :: jk, jl 
    371197      REAL(wp) :: zfice, ziter 
    372198      REAL(wp), DIMENSION(jpi,jpj) ::   z2d 
     
    379205      IF(lwp) THEN 
    380206         WRITE(numout,*) 
    381          WRITE(numout,*) 'ice_rst_read : read ice NetCDF restart file' 
    382          WRITE(numout,*) '~~~~~~~~~~~~~' 
     207         WRITE(numout,*) 'ice_rst_read: read ice NetCDF restart file' 
     208         WRITE(numout,*) '~~~~~~~~~~~~' 
    383209      ENDIF 
    384210 
     
    390216      IF(lwp) WRITE(numout,*) '   in any case we force it to nit000 - 1 : ', nit000 - 1 
    391217 
    392       !Control of date 
    393  
     218      ! Control of date 
    394219      IF( ( nit000 - NINT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 )   & 
    395220         &     CALL ctl_stop( 'ice_rst_read ===>>>> : problem with nit000 in ice restart',  & 
     
    451276      END DO 
    452277 
    453       CALL iom_get( numrir, jpdom_autoglo, 'u_ice'     , u_ice      ) 
    454       CALL iom_get( numrir, jpdom_autoglo, 'v_ice'     , v_ice      ) 
    455       CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  ) 
    456       CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  ) 
    457       CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i ) 
    458       CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass ) 
    459       CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) 
    460  
    461       ! In case Prather scheme is used for advection, read second order moments 
    462       ! ------------------------------------------------------------------------ 
    463       IF( ln_adv_Pra ) THEN 
    464  
    465          DO jl = 1, jpl  
    466             WRITE(zchar,'(I2.2)') jl 
    467             znam = 'sxice'//'_htc'//zchar 
    468             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    469             sxice(:,:,jl) = z2d(:,:) 
    470             znam = 'syice'//'_htc'//zchar 
    471             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    472             syice(:,:,jl) = z2d(:,:) 
    473             znam = 'sxxice'//'_htc'//zchar 
    474             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    475             sxxice(:,:,jl) = z2d(:,:) 
    476             znam = 'syyice'//'_htc'//zchar 
    477             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    478             syyice(:,:,jl) = z2d(:,:) 
    479             znam = 'sxyice'//'_htc'//zchar 
    480             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    481             sxyice(:,:,jl) = z2d(:,:) 
    482             znam = 'sxsn'//'_htc'//zchar 
    483             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    484             sxsn(:,:,jl) = z2d(:,:) 
    485             znam = 'sysn'//'_htc'//zchar 
    486             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    487             sysn(:,:,jl) = z2d(:,:) 
    488             znam = 'sxxsn'//'_htc'//zchar 
    489             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    490             sxxsn(:,:,jl) = z2d(:,:) 
    491             znam = 'syysn'//'_htc'//zchar 
    492             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    493             syysn(:,:,jl) = z2d(:,:) 
    494             znam = 'sxysn'//'_htc'//zchar 
    495             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    496             sxysn(:,:,jl) = z2d(:,:) 
    497             znam = 'sxa'//'_htc'//zchar 
    498             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    499             sxa(:,:,jl) = z2d(:,:) 
    500             znam = 'sya'//'_htc'//zchar 
    501             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    502             sya(:,:,jl) = z2d(:,:) 
    503             znam = 'sxxa'//'_htc'//zchar 
    504             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    505             sxxa(:,:,jl) = z2d(:,:) 
    506             znam = 'syya'//'_htc'//zchar 
    507             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    508             syya(:,:,jl) = z2d(:,:) 
    509             znam = 'sxya'//'_htc'//zchar 
    510             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    511             sxya(:,:,jl) = z2d(:,:) 
    512             znam = 'sxc0'//'_htc'//zchar 
    513             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    514             sxc0(:,:,jl) = z2d(:,:) 
    515             znam = 'syc0'//'_htc'//zchar 
    516             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    517             syc0(:,:,jl) = z2d(:,:) 
    518             znam = 'sxxc0'//'_htc'//zchar 
    519             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    520             sxxc0(:,:,jl) = z2d(:,:) 
    521             znam = 'syyc0'//'_htc'//zchar 
    522             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    523             syyc0(:,:,jl) = z2d(:,:) 
    524             znam = 'sxyc0'//'_htc'//zchar 
    525             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    526             sxyc0(:,:,jl) = z2d(:,:) 
    527             znam = 'sxsal'//'_htc'//zchar 
    528             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    529             sxsal(:,:,jl) = z2d(:,:) 
    530             znam = 'sysal'//'_htc'//zchar 
    531             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    532             sysal(:,:,jl) = z2d(:,:) 
    533             znam = 'sxxsal'//'_htc'//zchar 
    534             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    535             sxxsal(:,:,jl) = z2d(:,:) 
    536             znam = 'syysal'//'_htc'//zchar 
    537             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    538             syysal(:,:,jl) = z2d(:,:) 
    539             znam = 'sxysal'//'_htc'//zchar 
    540             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    541             sxysal(:,:,jl) = z2d(:,:) 
    542             znam = 'sxage'//'_htc'//zchar 
    543             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    544             sxage(:,:,jl) = z2d(:,:) 
    545             znam = 'syage'//'_htc'//zchar 
    546             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    547             syage(:,:,jl) = z2d(:,:) 
    548             znam = 'sxxage'//'_htc'//zchar 
    549             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    550             sxxage(:,:,jl) = z2d(:,:) 
    551             znam = 'syyage'//'_htc'//zchar 
    552             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    553             syyage(:,:,jl) = z2d(:,:) 
    554             znam = 'sxyage'//'_htc'//zchar 
    555             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    556             sxyage(:,:,jl)= z2d(:,:) 
    557          END DO 
    558          ! MV MP 2016 
    559          IF ( nn_pnd_scheme > 0 ) THEN 
    560             DO jl = 1, jpl  
    561                WRITE(zchar,'(I2.2)') jl 
    562                znam = 'sxap'//'_htc'//zchar 
    563                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    564                sxap(:,:,jl) = z2d(:,:) 
    565                znam = 'syap'//'_htc'//zchar 
    566                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    567                syap(:,:,jl) = z2d(:,:) 
    568                znam = 'sxxap'//'_htc'//zchar 
    569                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    570                sxxap(:,:,jl) = z2d(:,:) 
    571                znam = 'syyap'//'_htc'//zchar 
    572                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    573                syyap(:,:,jl) = z2d(:,:) 
    574                znam = 'sxyap'//'_htc'//zchar 
    575                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    576                sxyap(:,:,jl) = z2d(:,:) 
    577     
    578                znam = 'sxvp'//'_htc'//zchar 
    579                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    580                sxvp(:,:,jl) = z2d(:,:) 
    581                znam = 'syvp'//'_htc'//zchar 
    582                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    583                syvp(:,:,jl) = z2d(:,:) 
    584                znam = 'sxxvp'//'_htc'//zchar 
    585                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    586                sxxvp(:,:,jl) = z2d(:,:) 
    587                znam = 'syyvp'//'_htc'//zchar 
    588                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    589                syyvp(:,:,jl) = z2d(:,:) 
    590                znam = 'sxyvp'//'_htc'//zchar 
    591                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    592                sxyvp(:,:,jl) = z2d(:,:) 
    593             END DO 
    594          ENDIF 
    595          ! END MV MP 2016 
    596  
    597          CALL iom_get( numrir, jpdom_autoglo, 'sxopw ' ,  sxopw  ) 
    598          CALL iom_get( numrir, jpdom_autoglo, 'syopw ' ,  syopw  ) 
    599          CALL iom_get( numrir, jpdom_autoglo, 'sxxopw' ,  sxxopw ) 
    600          CALL iom_get( numrir, jpdom_autoglo, 'syyopw' ,  syyopw ) 
    601          CALL iom_get( numrir, jpdom_autoglo, 'sxyopw' ,  sxyopw ) 
    602  
    603          DO jl = 1, jpl  
    604             WRITE(zchar,'(I2.2)') jl 
    605             DO jk = 1, nlay_i  
    606                WRITE(zchar1,'(I2.2)') jk 
    607                znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 
    608                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    609                sxe(:,:,jk,jl) = z2d(:,:) 
    610                znam = 'sye'//'_il'//zchar1//'_htc'//zchar 
    611                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    612                sye(:,:,jk,jl) = z2d(:,:) 
    613                znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 
    614                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    615                sxxe(:,:,jk,jl) = z2d(:,:) 
    616                znam = 'syye'//'_il'//zchar1//'_htc'//zchar 
    617                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    618                syye(:,:,jk,jl) = z2d(:,:) 
    619                znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 
    620                CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    621                sxye(:,:,jk,jl) = z2d(:,:) 
    622             END DO 
    623          END DO 
    624          ! 
    625       END IF 
    626        
    627       ! clem: I do not understand why the following IF is needed 
    628       !       I suspect something inconsistent in the main code with option nn_icesal=1 
    629       IF( nn_icesal == 1 ) THEN 
    630          DO jl = 1, jpl  
    631             sm_i(:,:,jl) = rn_icesal 
    632             DO jk = 1, nlay_i  
    633                s_i(:,:,jk,jl) = rn_icesal 
    634             END DO 
    635          END DO 
    636       ENDIF 
    637       ! 
    638       !CALL iom_close( numrir ) !clem: closed in icestp.F90 
    639       ! 
     278      CALL iom_get( numrir, jpdom_autoglo, 'u_ice', u_ice ) 
     279      CALL iom_get( numrir, jpdom_autoglo, 'v_ice', v_ice ) 
     280 
    640281   END SUBROUTINE ice_rst_read 
    641282 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

    r8517 r8518  
    206206                                    CALL ice_wri( 1 )         ! -- Ice outputs  
    207207         ! 
    208          IF( kt == nit000 .AND. ln_rstart )   &                !!gm  I don't understand the ln_rstart, if needed, add a comment, please 
    209             &                       CALL iom_close( numrir )  ! close input ice restart file 
    210208         ! 
    211209         IF( lrst_ice )             CALL ice_rst_write( kt )  ! -- Ice restart file  
     
    251249      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing  
    252250      ierr = ierr + ice1D_alloc      ()      ! thermodynamics 
    253       ierr = ierr + ice_rdgrft_alloc ()      ! ridging/rafting 
    254251      ! 
    255252      IF( lk_mpp    )   CALL mpp_sum( ierr ) 
     
    258255      CALL ice_itd_init                ! ice thickness distribution initialization 
    259256      ! 
    260       IF( ln_icedyn )   THEN 
    261          CALL ice_dyn_init             ! set ice dynamics parameters 
    262          CALL ice_rdgrft_init          ! set ice ridging/rafting parameters 
    263          CALL ice_rhg_init             ! set ice rheology parameters 
    264          CALL ice_adv_init             ! set ice advection parameters 
    265       ENDIF 
    266  
    267       IF( ln_icethd ) THEN 
    268          CALL ice_thd_init             ! set ice thermodynics parameters 
    269          CALL ice_thd_sal_init         ! set ice salinity parameters 
    270       ENDIF 
    271     
    272257      ! MV MP 2016 
    273       CALL lim_mp_init                 ! set melt ponds parameters 
     258      CALL lim_mp_init                 ! set melt ponds parameters (clem: important to be located here) 
    274259      ! END MV MP 2016 
    275260      !                                ! Initial sea-ice state 
     
    283268      CALL ice_var_glo2eqv 
    284269      ! 
     270      IF( ln_icedyn ) THEN 
     271         CALL ice_dyn_init             ! set ice dynamics parameters 
     272         CALL ice_rdgrft_init          ! set ice ridging/rafting parameters 
     273         CALL ice_rhg_init             ! set ice rheology parameters 
     274         CALL ice_adv_init             ! set ice advection parameters 
     275      ENDIF 
     276 
     277      IF( ln_icethd ) THEN 
     278         CALL ice_thd_init             ! set ice thermodynics parameters 
     279         CALL ice_thd_sal_init         ! set ice salinity parameters 
     280      ENDIF    
     281      ! 
    285282      CALL ice_update_init             ! ice surface boundary condition 
    286283      ! 
     
    296293      ELSEWHERE                     ;   rn_amax_2d(:,:) = rn_amax_s  ! SH 
    297294      END WHERE 
     295 
     296      IF( ln_rstart )   CALL iom_close( numrir )  ! close input ice restart file 
    298297      ! 
    299298   END SUBROUTINE ice_init 
     
    429428      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
    430429      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
    431       hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
     430      hfx_err_rem(:,:) = 0._wp 
    432431      hfx_err_dif(:,:) = 0._wp 
    433432      wfx_err_sub(:,:) = 0._wp 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

    r8517 r8518  
    413413         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw          ) 
    414414         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub          ) 
    415          CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_err_1d (1:nidx), hfx_err          ) 
    416415         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res          ) 
    417416         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif   ) 
     
    499498         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw        ) 
    500499         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub        ) 
    501          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_1d (1:nidx), hfx_err        ) 
    502500         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res        ) 
    503501         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif ) 
     
    585583      ENDIF 
    586584      ! 
    587       IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'STOP', 'ice_thd_init : rn_hinew should be >= rn_himin' ) 
     585      IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'ice_thd_init : rn_hinew should be >= rn_himin' ) 
    588586      ! 
    589587      IF(lwp) WRITE(numout,*) 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90

    r8514 r8518  
    9494      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C  
    9595      REAL(wp) ::   ztmelt_i                  ! ice melting temperature 
    96       REAL(wp) ::   zhsu 
     96      REAL(wp) ::   z1_hsu 
    9797      REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
    9898      REAL(wp) ::   zdti_bnd = 1.e-4_wp       ! maximal authorized error on temperature  
     99      REAL(wp) ::   zcpi                      ! Ice specific heat 
     100      REAL(wp) ::   zhi                       !  
     101      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat 
    99102       
    100103      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow 
    101       REAL(wp), DIMENSION(jpij)     ::   ztsub       ! old surface temperature (before the iterative procedure ) 
    102       REAL(wp), DIMENSION(jpij)     ::   ztsubit     ! surface temperature at previous iteration 
     104      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration 
    103105      REAL(wp), DIMENSION(jpij)     ::   zh_i        ! ice layer thickness 
    104106      REAL(wp), DIMENSION(jpij)     ::   zh_s        ! snow layer thickness 
     
    106108      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface 
    107109      REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function 
    108       REAL(wp), DIMENSION(jpij)     ::   dzf        ! derivative of the surface flux function 
     110      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function 
    109111      REAL(wp), DIMENSION(jpij)     ::   zdti        ! current error on temperature 
    110       REAL(wp), DIMENSION(jpij)     ::   zdifcase    ! case of the equation resolution (1->4) 
    111112      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice 
    112       REAL(wp), DIMENSION(jpij)     ::   zihic 
    113113       
     114      REAL(wp), DIMENSION(jpij,nlay_i)   ::   z_i         ! Vertical cotes of the layers in the ice 
     115      REAL(wp), DIMENSION(jpij,nlay_s)   ::   z_s         ! Vertical cotes of the layers in the snow 
     116      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztib        ! Old temperature in the ice 
     117      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsb        ! Temporary temperature in the snow 
     118      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
     119      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
    114120      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity 
    115121      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
    116122      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice 
    117123      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice 
    118       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztib        ! Old temperature in the ice 
    119124      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice 
    120       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
    121       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zspeche_i   ! Ice specific heat 
    122       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   z_i         ! Vertical cotes of the layers in the ice 
    123125      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow 
    124126      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow 
    125127      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow 
    126128      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow 
    127       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
    128       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   ztsb        ! Temporary temperature in the snow 
    129       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   z_s         ! Vertical cotes of the layers in the snow 
    130129      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term 
    131130      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term 
    132131      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term 
    133132      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms 
    134       REAL(wp), DIMENSION(jpij)            ::   zdq, zq_ini, zhfx_err ! diag errors on heat 
     133      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat 
    135134      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
    136135       
     
    148147 
    149148      ! --- diag error on heat diffusion - PART 1 --- ! 
    150       zdq(:) = 0._wp ; zq_ini(:) = 0._wp       
    151149      DO ji = 1, nidx 
    152150         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & 
     
    167165      ! Ice / snow layers 
    168166      !-------------------- 
    169  
    170       z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer 
    171       z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer 
    172  
    173       DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer 
    174          DO ji = 1 , nidx 
    175             z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s 
    176          END DO 
    177       END DO 
    178  
    179       DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer 
    180          DO ji = 1 , nidx 
    181             z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i 
     167      z_s(1:nidx,1) = zh_s(1:nidx) 
     168      z_i(1:nidx,1) = zh_i(1:nidx) 
     169      DO jk = 2, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer 
     170         DO ji = 1 , nidx 
     171            z_s(ji,jk) = z_s(ji,jk-1) + zh_s(ji) 
     172         END DO 
     173      END DO 
     174 
     175      DO jk = 2, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer 
     176         DO ji = 1 , nidx 
     177            z_i(ji,jk) = z_i(ji,jk-1) + zh_i(ji) 
    182178         END DO 
    183179      END DO 
     
    187183      !------------------------------------------------------------------------------| 
    188184      ! 
    189       !------------------- 
    190       ! Computation of i0 
    191       !------------------- 
    192       ! i0 describes the fraction of solar radiation which does not contribute 
    193       ! to the surface energy budget but rather penetrates inside the ice. 
    194       ! We assume that no radiation is transmitted through the snow 
    195       ! If there is no no snow 
    196       ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface  
    197       ! zftrice = io.qsr_ice       is below the surface  
    198       ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
    199       ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 
    200       zhsu = 0.1_wp ! threshold for the computation of i0 
     185      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 
    201186      DO ji = 1 , nidx 
    202          ! switches 
    203          isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )  
    204          ! hs > 0, isnow = 1 
    205          zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )      
    206  
    207          i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 
    208       END DO 
    209  
    210       !------------------------------------------------------- 
    211       ! Solar radiation absorbed / transmitted at the surface 
    212       ! Derivative of the non solar flux 
    213       !------------------------------------------------------- 
    214       DO ji = 1 , nidx 
    215          zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface 
    216          zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer 
    217          dzf    (ji)    = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux  
    218          zqns_ice_b(ji) = qns_ice_1d(ji)                     ! store previous qns_ice_1d value 
     187         !------------------- 
     188         ! Computation of i0 
     189         !------------------- 
     190         ! i0 describes the fraction of solar radiation which does not contribute 
     191         ! to the surface energy budget but rather penetrates inside the ice. 
     192         ! We assume that no radiation is transmitted through the snow 
     193         ! If there is no no snow 
     194         ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface  
     195         ! zftrice = io.qsr_ice       is below the surface  
     196         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
     197         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 
     198         zhi = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )      
     199         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zhi * fr2_i0_1d(ji) ) 
     200 
     201         !------------------------------------------------------- 
     202         ! Solar radiation absorbed / transmitted at the surface 
     203         ! Derivative of the non solar flux 
     204         !------------------------------------------------------- 
     205         zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface 
     206         zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer 
     207         zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux  
     208         zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value 
    219209      END DO 
    220210 
     
    222212      ! Transmission - absorption of solar radiation in the ice 
    223213      !--------------------------------------------------------- 
    224  
    225       DO ji = 1, nidx           ! snow initialization 
    226          zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow 
    227       END DO 
    228  
    229       DO jk = 1, nlay_s          ! Radiation through snow 
     214      zradtr_s(1:nidx,0) = zftrice(1:nidx) 
     215      DO jk = 1, nlay_s 
    230216         DO ji = 1, nidx 
    231217            !                             ! radiation transmitted below the layer-th snow layer 
    232             zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) ) 
     218            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * z_s(ji,jk) ) 
    233219            !                             ! radiation absorbed by the layer-th snow layer 
    234220            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 
    235221         END DO 
    236222      END DO 
    237  
    238       DO ji = 1, nidx           ! ice initialization 
    239          zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 
    240       END DO 
    241  
    242       DO jk = 1, nlay_i          ! Radiation through ice 
     223          
     224      zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) ) 
     225      DO jk = 1, nlay_i  
    243226         DO ji = 1, nidx 
    244227            !                             ! radiation transmitted below the layer-th ice layer 
    245             zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 
     228            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * z_i(ji,jk) ) 
    246229            !                             ! radiation absorbed by the layer-th ice layer 
    247230            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 
     
    249232      END DO 
    250233 
    251       DO ji = 1, nidx           ! Radiation transmitted below the ice 
    252          ftr_ice_1d(ji) = zradtr_i(ji,nlay_i)  
    253       END DO 
     234      ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i)   ! record radiation transmitted below the ice 
    254235 
    255236      !------------------------------------------------------------------------------| 
     
    257238      !------------------------------------------------------------------------------| 
    258239      ! 
    259       DO ji = 1, nidx        ! Old surface temperature 
    260          ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr. 
    261          ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter 
    262          t_su_1d(ji) =  MIN( t_su_1d(ji), rt0 - ztsu_err )       ! necessary 
    263          zdti   (ji) =  1000._wp                                 ! initial value of error 
    264       END DO 
    265  
    266       DO jk = 1, nlay_s       ! Old snow temperature 
    267          DO ji = 1 , nidx 
    268             ztsb(ji,jk) =  t_s_1d(ji,jk) 
    269          END DO 
    270       END DO 
    271  
    272       DO jk = 1, nlay_i       ! Old ice temperature 
    273          DO ji = 1 , nidx 
    274             ztib(ji,jk) =  t_i_1d(ji,jk) 
    275          END DO 
    276       END DO 
     240      ztsub  (1:nidx)   =  t_su_1d(1:nidx)                          ! temperature at the previous iter 
     241      t_su_1d(1:nidx)   =  MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary 
     242      zdti   (1:nidx)   =  1000._wp                                 ! initial value of error 
     243      ztsb   (1:nidx,:) =  t_s_1d(1:nidx,:)                         ! Old snow temperature 
     244      ztib   (1:nidx,:) =  t_i_1d(1:nidx,:)                         ! Old ice temperature 
    277245 
    278246      iconv    =  0           ! number of iterations 
     
    289257         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T 
    290258            DO ji = 1 , nidx 
    291                ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
    292                ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
     259               ztcond_i(ji,0) = MAX( zkimin, rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) ) 
    293260            END DO 
    294261            DO jk = 1, nlay_i-1 
    295262               DO ji = 1 , nidx 
    296                   ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  & 
    297                      MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) 
    298                   ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 
     263                  ztcond_i(ji,jk) = MAX( zkimin, rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  & 
     264                     &                           MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) ) 
    299265               END DO 
    300266            END DO 
     
    302268         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 
    303269            DO ji = 1 , nidx 
    304                ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   & 
    305                   &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 )   
    306                ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
     270               ztcond_i(ji,0) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   & 
     271                  &                               - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) ) 
    307272            END DO 
    308273            DO jk = 1, nlay_i-1 
    309274               DO ji = 1 , nidx 
    310                   ztcond_i(ji,jk) = rcdic +                                                                       &  
     275                  ztcond_i(ji,jk) = MAX( zkimin, rcdic +                                                           &  
    311276                     &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              & 
    312277                     &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   & 
    313                      &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )   
    314                   ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 
     278                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )  ) 
    315279               END DO 
    316280            END DO 
    317281            DO ji = 1 , nidx 
    318                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   & 
    319                   &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )   
    320                ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 
     282               ztcond_i(ji,nlay_i) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   & 
     283                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) )   
    321284            END DO 
    322285         ENDIF 
     
    335298         SELECT CASE ( nn_monocat ) 
    336299 
    337          CASE (1,3) ! LIM3 
     300         CASE (1,3) 
    338301 
    339302            zepsilon = 0.1_wp 
     
    404367         DO jk = 1, nlay_i 
    405368            DO ji = 1 , nidx 
    406                ztitemp(ji,jk)   = t_i_1d(ji,jk) 
    407                zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 
    408                zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 ) 
     369               ztitemp(ji,jk) = t_i_1d(ji,jk) 
     370               zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 
     371               zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zcpi * zh_i(ji), epsi10 ) 
    409372            END DO 
    410373         END DO 
     
    425388            DO ji = 1 , nidx 
    426389               ! update of the non solar flux according to the update in T_su 
    427                qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 
     390               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
    428391            END DO 
    429392         ENDIF 
     
    507470                  !  case 1 : no surface melting - snow present                                  | 
    508471                  !------------------------------------------------------------------------------| 
    509                   zdifcase(ji)    =  1.0 
    510472                  numeqmin(ji)    =  1 
    511473                  numeqmax(ji)    =  nlay_i + nlay_s + 1 
     
    513475                  !!surface equation 
    514476                  ztrid(ji,1,1)  = 0.0 
    515                   ztrid(ji,1,2)  = dzf(ji) - zg1s * zkappa_s(ji,0) 
     477                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 
    516478                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0) 
    517                   zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji) 
     479                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji) 
    518480 
    519481                  !!first layer of snow equation 
     
    529491                  !------------------------------------------------------------------------------| 
    530492                  ! 
    531                   zdifcase(ji)    =  2.0 
    532493                  numeqmin(ji)    =  2 
    533494                  numeqmax(ji)    =  nlay_i + nlay_s + 1 
     
    552513                  !------------------------------------------------------------------------------| 
    553514                  ! 
    554                   zdifcase(ji)      =  3.0 
    555515                  numeqmin(ji)      =  nlay_s + 1 
    556516                  numeqmax(ji)      =  nlay_i + nlay_s + 1 
     
    558518                  !!surface equation    
    559519                  ztrid(ji,numeqmin(ji),1)   =  0.0 
    560                   ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1     
     520                  ztrid(ji,numeqmin(ji),2)   =  zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1     
    561521                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1 
    562                   zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji) 
     522                  zindterm(ji,numeqmin(ji))  =  zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji) 
    563523 
    564524                  !!first layer of ice equation 
     
    572532                  IF ( nlay_i == 1 ) THEN 
    573533                     ztrid(ji,numeqmin(ji),1)    =  0.0 
    574                      ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0 
     534                     ztrid(ji,numeqmin(ji),2)    =  zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 
    575535                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0 
    576536                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
     
    589549                  !------------------------------------------------------------------------------| 
    590550                  ! 
    591                   zdifcase(ji)    =  4.0 
    592551                  numeqmin(ji)    =  nlay_s + 2 
    593552                  numeqmax(ji)    =  nlay_i + nlay_s + 1 
     
    662621            ! surface temperature 
    663622            isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) ) 
    664             ztsubit(ji) = t_su_1d(ji) 
     623            ztsub(ji) = t_su_1d(ji) 
    665624            IF( t_su_1d(ji) < rt0 ) & 
    666625               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  & 
     
    676635         DO ji = 1 , nidx 
    677636            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  ) 
    678             zdti   (ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )      
     637            zdti   (ji) =  ABS( t_su_1d(ji) - ztsub(ji) )      
    679638         END DO 
    680639 
     
    754713 
    755714      ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 
     715      !     hfx_dif = Heat flux used to warm/cool ice in W.m-2 
     716      !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
    756717      DO ji = 1, nidx 
    757          zdq(ji)        = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & 
    758             &                              SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
     718         zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & 
     719            &                   SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
     720 
    759721         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
    760             zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice  
     722            zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)  
    761723         ELSE                          ! case T_su = 0degC 
    762             zhfx_err(ji) = fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice  
     724            zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
    763725         ENDIF 
    764          hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
     726         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 
    765727 
    766728         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 
    767          hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
    768       END DO  
    769  
    770       !----------------------------------------- 
    771       ! Heat flux used to warm/cool ice in W.m-2 
    772       !----------------------------------------- 
    773       DO ji = 1, nidx 
    774          IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
    775             hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   & 
    776                &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 
    777          ELSE                          ! case T_su = 0degC 
    778             hfx_dif_1d(ji) = hfx_dif_1d(ji) +    & 
    779                &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 
    780          ENDIF 
    781          ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
    782          hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji) 
     729         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 
     730 
    783731      END DO    
    784732      ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90

    r8514 r8518  
    2929   PUBLIC   ice_thd_sal        ! called by icethd module 
    3030   PUBLIC   ice_thd_sal_init   ! called by ice_init 
     31    
     32   ! ** namelist (namsal) ** 
     33   LOGICAL  ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F) 
     34   REAL(wp) ::   rn_sal_gd        !    restoring salinity for gravity drainage [PSU] 
     35   REAL(wp) ::   rn_time_gd       !    restoring time constant for gravity drainage (= 20 days) [s] 
     36   REAL(wp) ::   rn_sal_fl        !    restoring salinity for flushing [PSU] 
     37   REAL(wp) ::   rn_time_fl       !    restoring time constant for gravity drainage (= 10 days) [s] 
    3138 
    3239   !!---------------------------------------------------------------------- 
     
    115122      !! 
    116123      NAMELIST/namice_sal/ ln_icedS , nn_icesal , rn_icesal, rn_sal_gd, rn_time_gd,   & 
    117          &                rn_sal_fl, rn_time_fl, rn_simax , rn_simin  
     124         &                 rn_sal_fl, rn_time_fl, rn_simax , rn_simin  
    118125      !!------------------------------------------------------------------- 
    119126      ! 
     
    143150      ENDIF 
    144151      ! 
     152      IF( ln_icedS .AND. nn_icesal == 1 ) THEN 
     153         ln_icedS = .FALSE. 
     154         CALL ctl_warn('ln_icedS is set to false since constant ice salinity is chosen (nn_icesal=1)') 
     155      ENDIF 
     156      ! 
    145157   END SUBROUTINE ice_thd_sal_init 
    146158 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90

    r8514 r8518  
    3636!!gm end 
    3737   USE sbccpl         ! Surface boundary condition: coupled interface 
    38    USE icealb         ! albedo parameters 
     38   USE icealb         ! sea-ice: albedo parameters 
    3939   USE traqsr         ! add penetration of solar flux in the calculation of heat budget 
    40    USE domvvl         ! Variable volume 
    41    USE icectl         ! ??? 
     40   USE icectl         ! sea-ice: control prints 
    4241   USE bdy_oce , ONLY : ln_bdy 
    4342   ! 
     
    5857   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2] 
    5958   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean velocity   [m/s] 
    60    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   soce_0  , sice_0     ! cst SSS and ice salinity (levitating sea-ice)  
    6159 
    6260   !! * Substitutions 
     
    7371      !!             ***  ROUTINE ice_update_alloc *** 
    7472      !!------------------------------------------------------------------- 
    75       ALLOCATE( soce_0(jpi,jpj) , utau_oce(jpi,jpj) ,                       & 
    76          &      sice_0(jpi,jpj) , vtau_oce(jpi,jpj) , tmod_io(jpi,jpj), STAT=ice_update_alloc) 
    77          ! 
     73      ALLOCATE( utau_oce(jpi,jpj), vtau_oce(jpi,jpj), tmod_io(jpi,jpj), STAT=ice_update_alloc ) 
     74      ! 
    7875      IF( lk_mpp                )   CALL mpp_sum( ice_update_alloc ) 
    7976      IF( ice_update_alloc /= 0 )   CALL ctl_warn('ice_update_alloc: failed to allocate arrays') 
     77      ! 
    8078   END FUNCTION ice_update_alloc 
    8179 
     
    206204      ! 
    207205      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    208  
    209       !                    ! conservation test 
    210       IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final( 'iceupdate' ) 
    211       !                    ! control prints 
    212       IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 3, ' - Final state ice_update - ' ) 
    213       IF( ln_ctl    )   CALL ice_prt3D( 'iceupdate' ) 
    214       ! 
    215       IF( nn_timing == 1 )   CALL timing_stop('ice_update_flx') 
     206      ! 
     207      IF( lrst_ice ) THEN                       !* write snwice_mass fields in the restart file 
     208         CALL update_rst( 'WRITE', kt ) 
     209      ENDIF 
     210      ! 
     211      ! controls 
     212      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation 
     213      IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 
     214      IF( ln_ctl                         )   CALL ice_prt3D     ('iceupdate')                                       ! prints 
     215      IF( nn_timing == 1                 )   CALL timing_stop   ('ice_update_flx')                                  ! timing 
    216216      ! 
    217217   END SUBROUTINE ice_update_flx 
     
    320320      IF( ice_update_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_update_init : unable to allocate standard arrays' ) 
    321321      ! 
    322       soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating case 0 (i.e. virtual salt flux) 
    323       sice_0(:,:) = sice 
    324       WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &   ! reduced values in the Baltic Sea area 
    325          &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         )  
    326          soce_0(:,:) = 4._wp 
    327          sice_0(:,:) = 2._wp 
    328       END WHERE 
    329       ! 
    330       IF( .NOT.ln_rstart ) THEN              ! set  
    331          ! 
    332          snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )   ! snow+ice mass 
    333          snwice_mass_b(:,:) = snwice_mass(:,:) 
    334          ! 
    335          IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area 
    336             sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
    337             sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 
    338  
    339 !!gm I really don't like this stuff here...  Find a way to put that elsewhere or differently 
    340 !!gm 
    341             IF( .NOT.ln_linssh ) THEN 
    342                DO jk = 1,jpkm1                     ! adjust initial vertical scale factors 
    343                   e3t_n(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) )  ) 
    344                   e3t_b(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshb(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) )  ) 
    345                END DO 
    346                e3t_a(:,:,:) = e3t_b(:,:,:) 
    347 !!gm  we are in no-restart case, so sshn=sshb   ==>> faster calculation: 
    348 !!    REAL(wp) ::   ze3t   ! local scalar 
    349 !!    REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! workspace 
    350 !! 
    351 !!             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:) 
    352 !!             ELSEWHERE                ;   z2d(:,:) = 1._wp 
    353 !!             END WHERE 
    354 !!             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors                 
    355 !!                ze3t = e3t_0(:,:,jk) * z2d(:,:) 
    356 !!                e3t_n(:,:,jk) = ze3t 
    357 !!                e3t_b(:,:,jk) = ze3t 
    358 !!                e3t_a(:,:,jk) = ze3t 
    359 !!             END DO 
    360 !!gm  but since it is only done at the initialisation....  just the following can be acceptable: 
    361 !               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors 
    362 !                  e3t_n(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1))  ) 
    363 !               END DO 
    364 !               e3t_b(:,:,:) = e3t_n(:,:,:) 
    365 !               e3t_a(:,:,:) = e3t_n(:,:,:) 
    366 !!gm end                
    367                ! Reconstruction of all vertical scale factors at now and before time-steps 
    368                ! ========================================================================= 
    369                ! Horizontal scale factor interpolations 
    370                ! -------------------------------------- 
    371                CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 
    372                CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 
    373                CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    374                CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
    375                CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 
    376                ! Vertical scale factor interpolations 
    377                  ! ------------------------------------ 
    378                CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  ) 
    379                CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    380                CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
    381                CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    382                CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
    383                ! t- and w- points depth 
    384                ! ---------------------- 
    385 !!gm not sure of that.... 
    386                gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    387                gdepw_n(:,:,1) = 0.0_wp 
    388                gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
    389                DO jk = 2, jpk 
    390                   gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  ) 
    391                   gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) 
    392                   gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:) 
    393                END DO 
     322      CALL update_rst( 'READ' )  !* read or initialize all required files 
     323      ! 
     324   END SUBROUTINE ice_update_init 
     325 
     326   SUBROUTINE update_rst( cdrw, kt ) 
     327      !!--------------------------------------------------------------------- 
     328      !!                   ***  ROUTINE rhg_evp_rst  *** 
     329      !!                      
     330      !! ** Purpose :   Read or write RHG file in restart file 
     331      !! 
     332      !! ** Method  :   use of IOM library 
     333      !!---------------------------------------------------------------------- 
     334      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     335      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step 
     336      ! 
     337      INTEGER  ::   iter   ! local integer 
     338      INTEGER  ::   id1    ! local integer 
     339      !!---------------------------------------------------------------------- 
     340      ! 
     341      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize 
     342         !                                   ! --------------- 
     343         IF( ln_rstart ) THEN                   !* Read the restart file 
     344            ! 
     345            id1 = iom_varid( numrir, 'snwice_mass' , ldstop = .FALSE. ) 
     346            ! 
     347            IF( id1 > 0 ) THEN                       ! fields exist 
     348               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass   ) 
     349               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) 
     350            ELSE                                     ! start from rest 
     351               IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it' 
     352               snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 
     353               snwice_mass_b(:,:) = snwice_mass(:,:) 
    394354            ENDIF 
     355         ELSE                                   !* Start from rest 
     356            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass' 
     357            snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 
     358            snwice_mass_b(:,:) = snwice_mass(:,:) 
    395359         ENDIF 
    396       ENDIF ! .NOT. ln_rstart 
    397       ! 
    398    END SUBROUTINE ice_update_init 
     360         ! 
     361      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
     362         !                                   ! ------------------- 
     363         IF(lwp) WRITE(numout,*) '---- update-rst ----' 
     364         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1 
     365         ! 
     366         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass   ) 
     367         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) 
     368         ! 
     369      ENDIF 
     370      ! 
     371   END SUBROUTINE update_rst 
    399372 
    400373#else 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90

    r8500 r8518  
    192192      CALL iom_put ('hfxsnw'     , hfx_snw(:,:)         )   !   
    193193      CALL iom_put ('hfxsub'     , hfx_sub(:,:)         )   !   
    194       CALL iom_put ('hfxerr'     , hfx_err(:,:)         )   !   
     194      CALL iom_put ('hfxerr'     , hfx_err_dif(:,:)     )   !   
    195195      CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)     )   !   
    196196       
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90

    r8226 r8518  
    8888   INTEGER ::   nitrst                !: time step at which restart file should be written 
    8989   LOGICAL ::   lrst_oce              !: logical to control the oce restart write  
     90   LOGICAL ::   lrst_ice              !: logical to control the ice restart write  
    9091   INTEGER ::   numror = 0            !: logical unit for ocean restart (read). Init to 0 is needed for SAS (in daymod.F90) 
     92   INTEGER ::   numrir                !: logical unit for ice   restart (read) 
    9193   INTEGER ::   numrow                !: logical unit for ocean restart (write) 
     94   INTEGER ::   numriw                !: logical unit for ice   restart (write) 
    9295   INTEGER ::   nrst_lst              !: number of restart to output next 
    9396 
Note: See TracChangeset for help on using the changeset viewer.