Changeset 11381


Ignore:
Timestamp:
2019-07-31T16:44:56+02:00 (14 months ago)
Author:
girrmann
Message:

dev_r10984_HPC-13 : repair faulty commit [11380], see #2308

Location:
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/cfgs/SHARED/namelist_ref

    r11365 r11381  
    13061306   jpni        =   0       !  jpni   number of processors following i (set automatically if < 1) 
    13071307   jpnj        =   0       !  jpnj   number of processors following j (set automatically if < 1) 
     1308   nn_hlts     =   1       !  added halo width for time splitting 
    13081309/ 
    13091310!----------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/ice.F90

    r11380 r11381  
    102102   !! vt_i        |      -      |    Total ice vol. per unit area | m     | 
    103103   !! vt_s        |      -      |    Total snow vol. per unit ar. | m     | 
     104   !! st_i        |      -      |    Total Sea ice salt content   | pss.m | 
    104105   !! sm_i        |      -      |    Mean sea ice salinity        | pss   | 
    105106   !! tm_i        |      -      |    Mean sea ice temperature     | K     | 
     
    135136   REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice 
    136137   LOGICAL , PUBLIC ::   ln_landfast_L16  !: landfast ice parameterizationfrom lemieux2016  
    137    LOGICAL , PUBLIC ::   ln_landfast_home !: landfast ice parameterizationfrom home made  
    138138   REAL(wp), PUBLIC ::   rn_depfra        !:    fraction of ocean depth that ice must reach to initiate landfast ice 
    139139   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home)  
     
    252252   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_err_sub !: mass flux error after sublimation                        [kg.m-2.s-1] 
    253253 
    254    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total)        [s-1] 
    255  
    256254   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice bottom growth                   [pss.kg.m-2.s-1 => g.m-2.s-1] 
    257255   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bom     !: salt flux due to ice bottom melt                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
     
    309307   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice !: components of the ice velocity                          (m/s) 
    310308   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i , vt_s  !: ice and snow total volume per unit area                 (m) 
     309   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   st_i         !: Total ice salinity content                              (pss.m) 
    311310   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i         !: ice total fractional area (ice concentration) 
    312311   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ato_i        !: =1-at_i ; total open water fractional area 
     
    409408         &      wfx_bog    (jpi,jpj) , wfx_dyn   (jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,           & 
    410409         &      wfx_res    (jpi,jpj) , wfx_sni   (jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,           & 
    411          &      afx_tot    (jpi,jpj) , rn_amax_2d(jpi,jpj),                                                  & 
     410         &      rn_amax_2d (jpi,jpj) ,                                                                       & 
    412411         &      qsb_ice_bot(jpi,jpj) , qlead     (jpi,jpj) ,                                                 & 
    413412         &      sfx_res    (jpi,jpj) , sfx_bri   (jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  & 
     
    429428      ii = ii + 1 
    430429      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                   & 
    431          &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) ,  & 
    432          &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s (jpi,jpj) ,  & 
    433          &      sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s (jpi,jpj) ,  & 
     430         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , st_i(jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) ,  & 
     431         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s(jpi,jpj) ,  & 
     432         &      sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s(jpi,jpj) ,  & 
    434433         &      om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj)            , STAT=ierr(ii) ) 
    435434 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icecor.F90

    r11380 r11381  
    8484      !                             !----------------------------------------------------- 
    8585      IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        ! 
    86       !                             !----------------------------------------------------- 
     86         !                          !----------------------------------------------------- 
    8787         zzc = rhoi * r1_rdtice 
    8888         DO jl = 1, jpl 
     
    117117            END DO 
    118118         END DO 
    119          CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1., v_ice, 'V', -1. )            ! lateral boundary conditions 
     119         CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1., v_ice, 'V', -1. ) 
    120120      ENDIF 
    121121 
    122 !!gm I guess the trends are only out on demand  
    123 !!   So please, only do this is it exite an iom_use of on a these variables 
    124 !!   furthermore, only allocate the diag_ arrays in this case  
    125 !!   and do the iom_put here so that it is only a local allocation 
    126 !!gm  
    127122      !                             !----------------------------------------------------- 
    128123      SELECT CASE( kn )             !  Diagnostics                                       ! 
     
    143138         END DO 
    144139         !                       ! concentration tendency (dynamics) 
    145          zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice  
    146          afx_tot(:,:) = zafx(:,:) 
    147          IF( iom_use('afxdyn') )   CALL iom_put( 'afxdyn' , zafx(:,:) ) 
     140         IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN  
     141            zafx(:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice  
     142            CALL iom_put( 'afxdyn' , zafx ) 
     143         ENDIF 
    148144         ! 
    149145      CASE( 2 )                        !--- thermo trend diagnostics & ice aging 
     
    164160         END DO 
    165161         !                       ! concentration tendency (total + thermo) 
    166          zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 
    167          afx_tot(:,:) = afx_tot(:,:) + zafx(:,:) 
    168          IF( iom_use('afxthd') )   CALL iom_put( 'afxthd' , zafx(:,:) ) 
    169          IF( iom_use('afxtot') )   CALL iom_put( 'afxtot' , afx_tot(:,:) ) 
     162         IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN  
     163            zafx(:,:) = zafx(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 
     164            CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice ) 
     165            CALL iom_put( 'afxtot' , zafx ) 
     166         ENDIF 
    170167         ! 
    171168      END SELECT 
     
    174171      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    175172      IF( ln_ctl         )   CALL ice_prt3D   ('icecor')                                                             ! prints 
    176       IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                   ! prints 
     173      IF( ln_icectl .AND. kn == 2 ) & 
     174         &                   CALL ice_prt     ( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                       ! prints 
    177175      IF( ln_timing      )   CALL timing_stop ('icecor')                                                             ! timing 
    178176      ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedia.F90

    r11380 r11381  
    3434   PUBLIC   ice_dia_init   ! called in icestp.F90 
    3535 
    36    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents 
     36   REAL(wp), SAVE ::   z1_e1e2  ! inverse of the ocean area 
     37   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini                    ! initial volume, salt and heat contents 
    3738   REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends 
    3839    
     
    8081      ENDIF 
    8182 
    82 !!gm glob_sum includes a " * tmask_i ", so remove  " * tmask(:,:,1) " 
    83  
     83      IF( kt == nit000 ) THEN 
     84         z1_e1e2 = 1._wp / glob_sum( 'icedia', e1e2t(:,:) ) 
     85      ENDIF 
     86       
    8487      ! ----------------------- ! 
    85       ! 1 -  Contents ! 
     88      ! 1 -  Contents           ! 
    8689      ! ----------------------- ! 
    87       zbg_ivol = glob_sum( 'icedia', vt_i(:,:) * e1e2t(:,:) ) * 1.e-9                  ! ice volume (km3) 
    88       zbg_svol = glob_sum( 'icedia', vt_s(:,:) * e1e2t(:,:) ) * 1.e-9                  ! snow volume (km3) 
    89       zbg_area = glob_sum( 'icedia', at_i(:,:) * e1e2t(:,:) ) * 1.e-6                  ! area (km2) 
    90       zbg_isal = glob_sum( 'icedia', SUM( sv_i(:,:,:), dim=3 ) * e1e2t(:,:) ) * 1.e-9  ! salt content (pss*km3) 
    91       zbg_item = glob_sum( 'icedia', et_i * e1e2t(:,:) ) * 1.e-20                      ! heat content (1.e20 J) 
    92       zbg_stem = glob_sum( 'icedia', et_s * e1e2t(:,:) ) * 1.e-20                      ! heat content (1.e20 J) 
    93        
     90      IF(  iom_use('ibgvol_tot' ) .OR. iom_use('sbgvol_tot' ) .OR. iom_use('ibgarea_tot') .OR. & 
     91         & iom_use('ibgsalt_tot') .OR. iom_use('ibgheat_tot') .OR. iom_use('sbgheat_tot') ) THEN 
     92 
     93         zbg_ivol = glob_sum( 'icedia', vt_i(:,:) * e1e2t(:,:) ) * 1.e-9  ! ice volume (km3) 
     94         zbg_svol = glob_sum( 'icedia', vt_s(:,:) * e1e2t(:,:) ) * 1.e-9  ! snow volume (km3) 
     95         zbg_area = glob_sum( 'icedia', at_i(:,:) * e1e2t(:,:) ) * 1.e-6  ! area (km2) 
     96         zbg_isal = glob_sum( 'icedia', st_i(:,:) * e1e2t(:,:) ) * 1.e-9  ! salt content (pss*km3) 
     97         zbg_item = glob_sum( 'icedia', et_i(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J) 
     98         zbg_stem = glob_sum( 'icedia', et_s(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J) 
     99 
     100         CALL iom_put( 'ibgvol_tot'  , zbg_ivol )  
     101         CALL iom_put( 'sbgvol_tot'  , zbg_svol )  
     102         CALL iom_put( 'ibgarea_tot' , zbg_area )  
     103         CALL iom_put( 'ibgsalt_tot' , zbg_isal )  
     104         CALL iom_put( 'ibgheat_tot' , zbg_item )  
     105         CALL iom_put( 'sbgheat_tot' , zbg_stem )  
     106  
     107      ENDIF 
     108 
    94109      ! ---------------------------! 
    95110      ! 2 - Trends due to forcing  ! 
    96111      ! ---------------------------! 
     112      ! they must be kept outside an IF(iom_use) because of the call to dia_rst below 
    97113      z_frc_volbot = r1_rau0 * glob_sum( 'icedia', -( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-ocean  
    98114      z_frc_voltop = r1_rau0 * glob_sum( 'icedia', -( wfx_sub(:,:) + wfx_spr(:,:) )                    * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-atm 
     
    106122      frc_temtop  = frc_temtop  + z_frc_temtop  * rdt_ice ! 1.e20 J 
    107123      frc_tembot  = frc_tembot  + z_frc_tembot  * rdt_ice ! 1.e20 J 
     124 
     125      CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water)  
     126      CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water)  
     127      CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)    
     128      CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)    
     129      CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)    
     130 
     131      IF(  iom_use('ibgfrchfxtop') .OR. iom_use('ibgfrchfxbot') ) THEN 
     132         CALL iom_put( 'ibgfrchfxtop' , frc_temtop * z1_e1e2 * 1.e-20 * kt*rdt ) ! heat on top of ice/snw/ocean      (W/m2) 
     133         CALL iom_put( 'ibgfrchfxbot' , frc_tembot * z1_e1e2 * 1.e-20 * kt*rdt ) ! heat on top of ocean(below ice)   (W/m2)  
     134      ENDIF 
     135       
     136      ! ---------------------------------- ! 
     137      ! 3 -  Content variations and drifts ! 
     138      ! ---------------------------------- ! 
     139      IF(  iom_use('ibgvolume') .OR. iom_use('ibgsaltco') .OR. iom_use('ibgheatco') .OR. iom_use('ibgheatfx') ) THEN 
    108140             
    109       ! ----------------------- ! 
    110       ! 3 -  Content variations ! 
    111       ! ----------------------- ! 
    112       zdiff_vol = r1_rau0 * glob_sum( 'icedia', ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater trend (km3)  
    113       zdiff_sal = r1_rau0 * glob_sum( 'icedia', ( rhoi* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! salt content trend (km3*pss) 
    114       zdiff_tem =           glob_sum( 'icedia', ( et_i(:,:) + et_s(:,:)           - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20  ! heat content trend (1.e20 J) 
    115       !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check) 
    116  
    117       ! ----------------------- ! 
    118       ! 4 -  Drifts             ! 
    119       ! ----------------------- ! 
    120       zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot ) 
    121       zdiff_sal = zdiff_sal - frc_sal 
    122       zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop ) 
    123  
    124       ! ----------------------- ! 
    125       ! 5 - Diagnostics writing ! 
    126       ! ----------------------- ! 
    127 !!gm I don't understand the division by the ocean surface (i.e. glob_sum( 'icedia', e1e2t(:,:) ) * 1.e-20 * kt*rdt ) 
    128 !!   and its multiplication bu kt ! is it really what we want ? what is this quantity ? 
    129 !!   IF it is really what we want, compute it at kt=nit000, not 3 time by time-step ! 
    130 !!   kt*rdt  : you mean rdtice ? 
    131 !!gm 
    132       ! 
    133       IF( iom_use('ibgvolume')    )   CALL iom_put( 'ibgvolume' , zdiff_vol     )   ! ice/snow volume  drift            (km3 equivalent ocean water)          
    134       IF( iom_use('ibgsaltco')    )   CALL iom_put( 'ibgsaltco' , zdiff_sal     )   ! ice salt content drift            (psu*km3 equivalent ocean water) 
    135       IF( iom_use('ibgheatco')    )   CALL iom_put( 'ibgheatco' , zdiff_tem     )   ! ice/snow heat content drift       (1.e20 J) 
    136       IF( iom_use('ibgheatfx')    )   CALL iom_put( 'ibgheatfx' ,               &   ! ice/snow heat flux drift          (W/m2) 
    137          &                                                     zdiff_tem /glob_sum( 'icedia', e1e2t(:,:) * 1.e-20 * kt*rdt ) ) 
    138  
    139       IF( iom_use('ibgfrcvoltop') )   CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water)  
    140       IF( iom_use('ibgfrcvolbot') )   CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water)  
    141       IF( iom_use('ibgfrcsal')    )   CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)    
    142       IF( iom_use('ibgfrctemtop') )   CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)    
    143       IF( iom_use('ibgfrctembot') )   CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)    
    144       IF( iom_use('ibgfrchfxtop') )   CALL iom_put( 'ibgfrchfxtop' ,            &   ! heat on top of ice/snw/ocean      (W/m2)  
    145          &                                                          frc_temtop / glob_sum( 'icedia', e1e2t(:,:) ) * 1.e-20 * kt*rdt  ) 
    146       IF( iom_use('ibgfrchfxbot') )   CALL iom_put( 'ibgfrchfxbot' ,            &   ! heat on top of ocean(below ice)   (W/m2)  
    147          &                                                          frc_tembot / glob_sum( 'icedia', e1e2t(:,:) ) * 1.e-20 * kt*rdt  ) 
    148  
    149       IF( iom_use('ibgvol_tot' )  )   CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                       (km3) 
    150       IF( iom_use('sbgvol_tot' )  )   CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                      (km3) 
    151       IF( iom_use('ibgarea_tot')  )   CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                         (km2) 
    152       IF( iom_use('ibgsalt_tot')  )   CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content             (pss*km3) 
    153       IF( iom_use('ibgheat_tot')  )   CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                 (1.e20 J) 
    154       IF( iom_use('sbgheat_tot')  )   CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                (1.e20 J) 
    155       ! 
     141         zdiff_vol = r1_rau0 * glob_sum( 'icedia', ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater trend (km3)  
     142         zdiff_sal = r1_rau0 * glob_sum( 'icedia', ( rhoi*st_i(:,:)                  - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! salt content trend (km3*pss) 
     143         zdiff_tem =           glob_sum( 'icedia', ( et_i(:,:) + et_s(:,:)           - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20  ! heat content trend (1.e20 J) 
     144         !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check) 
     145          
     146         zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot ) 
     147         zdiff_sal = zdiff_sal - frc_sal 
     148         zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop ) 
     149          
     150         CALL iom_put( 'ibgvolume' , zdiff_vol )   ! ice/snow volume  drift            (km3 equivalent ocean water)          
     151         CALL iom_put( 'ibgsaltco' , zdiff_sal )   ! ice salt content drift            (psu*km3 equivalent ocean water) 
     152         CALL iom_put( 'ibgheatco' , zdiff_tem )   ! ice/snow heat content drift       (1.e20 J) 
     153         ! 
     154      ENDIF 
     155       
    156156      IF( lrst_ice )   CALL ice_dia_rst( 'WRITE', kt_ice ) 
    157157      ! 
     
    248248            vol_loc_ini(:,:) = rhoi * vt_i(:,:) + rhos * vt_s(:,:)  ! ice/snow volume (kg/m2) 
    249249            tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                ! ice/snow heat content (J) 
    250             sal_loc_ini(:,:) = rhoi * SUM( sv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2) 
     250            sal_loc_ini(:,:) = rhoi * st_i(:,:)                     ! ice salt content (pss*kg/m2) 
    251251         ENDIF 
    252252         ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedyn.F90

    r11380 r11381  
    163163            END DO 
    164164            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
    165             CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     165            ! output 
     166            CALL iom_put( 'icediv' , zdivu_i ) 
     167 
    166168            DEALLOCATE( zdivu_i ) 
    167169 
     
    219221      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  & 
    220222         &             rn_ishlat ,                                                           & 
    221          &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 
     223         &             ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 
    222224      !!------------------------------------------------------------------- 
    223225      ! 
     
    242244         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat 
    243245         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16 
    244          WRITE(numout,*) '      Landfast: param from home made                         ln_landfast_home= ', ln_landfast_home 
    245246         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra 
    246247         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr 
     
    269270      ENDIF 
    270271      !                                      !--- Landfast ice 
    271       IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp 
    272       ! 
    273       IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN 
    274          CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' ) 
    275       ENDIF 
     272      IF( .NOT.ln_landfast_L16 )   tau_icebfr(:,:) = 0._wp 
    276273      ! 
    277274      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedyn_adv.F90

    r11380 r11381  
    100100      diag_trp_vi(:,:) = SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_rdtice 
    101101      diag_trp_vs(:,:) = SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_rdtice 
    102       IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoi          )   ! ice mass transport 
    103       IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhos          )   ! snw mass transport 
    104       IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoi * 1.e-03 )   ! salt mass transport (kg/m2/s) 
    105       IF( iom_use('dihctrp') )   CALL iom_put( "dihctrp" , -diag_trp_ei                )   ! advected ice heat content (W/m2) 
    106       IF( iom_use('dshctrp') )   CALL iom_put( "dshctrp" , -diag_trp_es                )   ! advected snw heat content (W/m2) 
     102      IF( iom_use('icemtrp') )   CALL iom_put( 'icemtrp' , diag_trp_vi * rhoi          )   ! ice mass transport 
     103      IF( iom_use('snwmtrp') )   CALL iom_put( 'snwmtrp' , diag_trp_vs * rhos          )   ! snw mass transport 
     104      IF( iom_use('salmtrp') )   CALL iom_put( 'salmtrp' , diag_trp_sv * rhoi * 1.e-03 )   ! salt mass transport (kg/m2/s) 
     105      IF( iom_use('dihctrp') )   CALL iom_put( 'dihctrp' , -diag_trp_ei                 )   ! advected ice heat content (W/m2) 
     106      IF( iom_use('dshctrp') )   CALL iom_put( 'dshctrp' , -diag_trp_es                 )   ! advected snw heat content (W/m2) 
    107107 
    108108      ! controls 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedyn_rhg.F90

    r11380 r11381  
    6969         WRITE(numout,*)'ice_dyn_rhg: sea-ice rheology' 
    7070         WRITE(numout,*)'~~~~~~~~~~~' 
    71       ENDIF 
    72       ! 
    73       IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization 
    74          tau_icebfr(:,:) = 0._wp 
    75          DO jl = 1, jpl 
    76             WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
    77          END DO 
    7871      ENDIF 
    7972      ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedyn_rhg_evp.F90

    r11380 r11381  
    112112      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    113113      !! 
    114       LOGICAL, PARAMETER ::   ll_bdy_substep = .TRUE. ! temporary option to call bdy at each sub-time step (T) 
    115       !                                                                              or only at the main time step (F) 
    116114      INTEGER ::   ji, jj       ! dummy loop indices 
    117115      INTEGER ::   jter         ! local integers 
     
    137135      ! 
    138136      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points 
    139       REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points 
     137      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points 
    140138      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points 
    141139      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    142       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points 
    143       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ib , ztauV_ib             ! ice-bottom stress at U-V points (landfast param) 
    144       REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points 
    145140      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
    146       REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses 
    147141      ! 
    148142      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
     
    152146      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
    153147      !                                                                 !    ice bottom surface if ice is embedded    
    154       REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array 
    155       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array 
    156       ! 
    157       REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays 
    158       REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence 
     148      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses 
     149      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points 
     150      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array 
     151      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points 
     152      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points 
     153      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast) 
     154      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast) 
     155      ! 
     156      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays 
     157      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence 
    159158      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice 
    160159 
     
    163162      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
    164163      !! --- diags 
    165       REAL(wp), DIMENSION(jpi,jpj) ::   zswi 
     164      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00 
    166165      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3 
    167166      !! --- SIMIP diags 
    168       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice    
    169       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice 
    170       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2) 
    171       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2) 
    172       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2) 
    173       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2) 
    174       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2) 
    175       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2) 
    176       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress 
    177       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress   
    178167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
    179168      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) 
     
    264253 
    265254      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
    266       IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN   ;   zkt = rn_tensile 
    267       ELSE                                               ;   zkt = 0._wp 
     255      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile 
     256      ELSE                         ;   zkt = 0._wp 
    268257      ENDIF 
    269258      ! 
     
    308297             
    309298            ! Drag ice-atm. 
    310             zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
    311             zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
     299            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
     300            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
    312301 
    313302            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
     
    316305 
    317306            ! masks 
    318             zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
    319             zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
     307            zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
     308            zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
    320309 
    321310            ! switches 
    322             IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zswitchU(ji,jj) = 0._wp 
    323             ELSE                                                   ;   zswitchU(ji,jj) = 1._wp   ;   ENDIF 
    324             IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zswitchV(ji,jj) = 0._wp 
    325             ELSE                                                   ;   zswitchV(ji,jj) = 1._wp   ;   ENDIF 
     311            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
     312            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
     313            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
     314            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF 
    326315 
    327316         END DO 
     
    339328               ! ice-bottom stress at U points 
    340329               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 
    341                zTauU_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
     330               ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    342331               ! ice-bottom stress at V points 
    343332               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 
    344                zTauV_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
     333               ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    345334               ! ice_bottom stress at T points 
    346335               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 
    347                tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
     336               tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    348337            END DO 
    349338         END DO 
    350339         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 
    351340         ! 
    352       ELSEIF( ln_landfast_home ) THEN          !-- Home made 
     341      ELSE                               !-- no landfast 
    353342         DO jj = 2, jpjm1 
    354343            DO ji = fs_2, fs_jpim1 
    355                zTauU_ib(ji,jj) = tau_icebfr(ji,jj) 
    356                zTauV_ib(ji,jj) = tau_icebfr(ji,jj) 
    357             END DO 
    358          END DO 
    359          ! 
    360       ELSE                                     !-- no landfast 
    361          DO jj = 2, jpjm1 
    362             DO ji = fs_2, fs_jpim1 
    363                zTauU_ib(ji,jj) = 0._wp 
    364                zTauV_ib(ji,jj) = 0._wp 
     344               ztaux_base(ji,jj) = 0._wp 
     345               ztauy_base(ji,jj) = 0._wp 
    365346            END DO 
    366347         END DO 
    367348      ENDIF 
    368       IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) ) 
    369349 
    370350      !------------------------------------------------------------------------------! 
     
    504484                  !                 !--- tau_bottom/v_ice 
    505485                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
    506                   zTauB = - zTauV_ib(ji,jj) / zvel 
     486                  zTauB = ztauy_base(ji,jj) / zvel 
     487                  !                 !--- OceanBottom-to-Ice stress 
     488                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
    507489                  ! 
    508490                  !                 !--- Coriolis at V-points (energy conserving formulation) 
    509                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     491                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    510492                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    511493                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    512494                  ! 
    513495                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    514                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     496                  zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    515497                  ! 
    516498                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    517                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     499                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    518500                  ! 
    519501                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    520                   v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
    521                      &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    522                      &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    523                                     + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    524                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    525                      &           ) * zmaskV(ji,jj) 
     502                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
     503                        &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     504                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     505                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     506                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     507                        &           )   * zmsk00y(ji,jj) 
    526508                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    527                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             & ! previous velocity 
    528                      &                                     + zTauE + zTauO * v_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    529                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    530                      &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    531                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    532                      &            ) * zmaskV(ji,jj) 
     509                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
     510                        &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     511                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     512                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     513                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     514                        &            )   * zmsk00y(ji,jj) 
    533515                  ENDIF 
    534516               END DO 
     
    540522            CALL agrif_interp_ice( 'V' ) 
    541523#endif 
    542             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' ) 
     524            IF( ln_bdy CALL bdy_ice_dyn( 'V' ) 
    543525            ! 
    544526            DO jj = 2, jpjm1 
     
    552534                  !                 !--- tau_bottom/u_ice 
    553535                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    554                   zTauB = - zTauU_ib(ji,jj) / zvel 
     536                  zTauB = ztaux_base(ji,jj) / zvel 
     537                  !                 !--- OceanBottom-to-Ice stress 
     538                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 
    555539                  ! 
    556540                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    557                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
     541                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    558542                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    559543                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    560544                  ! 
    561545                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    562                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     546                  zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    563547                  ! 
    564548                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    565                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     549                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    566550                  ! 
    567551                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    568                   u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
    569                      &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    570                      &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    571                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
    572                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    573                      &            ) * zmaskU(ji,jj) 
     552                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
     553                        &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     554                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     555                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     556                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     557                        &            )   * zmsk00x(ji,jj) 
    574558                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    575                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             & ! previous velocity 
    576                      &                                     + zTauE + zTauO * u_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    577                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    578                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    579                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    580                      &            ) * zmaskU(ji,jj) 
     559                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
     560                        &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     561                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     562                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     563                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     564                        &            )   * zmsk00x(ji,jj) 
    581565                  ENDIF 
    582566               END DO 
     
    588572            CALL agrif_interp_ice( 'U' ) 
    589573#endif 
    590             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' ) 
     574            IF( ln_bdy CALL bdy_ice_dyn( 'U' ) 
    591575            ! 
    592576         ELSE ! odd iterations 
     
    602586                  !                 !--- tau_bottom/u_ice 
    603587                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    604                   zTauB = - zTauU_ib(ji,jj) / zvel 
     588                  zTauB = ztaux_base(ji,jj) / zvel 
     589                  !                 !--- OceanBottom-to-Ice stress 
     590                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 
    605591                  ! 
    606592                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    607                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
     593                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    608594                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    609595                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    610596                  ! 
    611597                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    612                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     598                  zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    613599                  ! 
    614600                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    615                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     601                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    616602                  ! 
    617603                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    618                   u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
    619                      &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    620                      &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    621                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
    622                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    623                      &            ) * zmaskU(ji,jj) 
     604                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
     605                        &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     606                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     607                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     608                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     609                        &            )   * zmsk00x(ji,jj) 
    624610                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    625                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             & ! previous velocity 
    626                      &                                     + zTauE + zTauO * u_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    627                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    628                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    629                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    630                      &            ) * zmaskU(ji,jj) 
     611                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
     612                        &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     613                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     614                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     615                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     616                        &            )   * zmsk00x(ji,jj) 
    631617                  ENDIF 
    632618               END DO 
     
    638624            CALL agrif_interp_ice( 'U' ) 
    639625#endif 
    640             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' ) 
     626            IF( ln_bdy CALL bdy_ice_dyn( 'U' ) 
    641627            ! 
    642628            DO jj = 2, jpjm1 
     
    650636                  !                 !--- tau_bottom/v_ice 
    651637                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
    652                   zTauB = - zTauV_ib(ji,jj) / zvel 
     638                  zTauB = ztauy_base(ji,jj) / zvel 
     639                  !                 !--- OceanBottom-to-Ice stress 
     640                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
    653641                  ! 
    654642                  !                 !--- Coriolis at v-points (energy conserving formulation) 
    655                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     643                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    656644                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    657645                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    658646                  ! 
    659647                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    660                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     648                  zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    661649                  ! 
    662650                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    663                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     651                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    664652                  ! 
    665653                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    666                   v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
    667                      &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    668                      &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    669                                     + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    670                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    671                      &           ) * zmaskV(ji,jj) 
     654                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
     655                        &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     656                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     657                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     658                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     659                        &           )   * zmsk00y(ji,jj) 
    672660                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    673                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             & ! previous velocity 
    674                      &                                     + zTauE + zTauO * v_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    675                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    676                      &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    677                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    678                      &            ) * zmaskV(ji,jj) 
     661                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
     662                        &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     663                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     664                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     665                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     666                        &            )   * zmsk00y(ji,jj) 
    679667                  ENDIF 
    680668               END DO 
     
    686674            CALL agrif_interp_ice( 'V' ) 
    687675#endif 
    688             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' ) 
     676            IF( ln_bdy CALL bdy_ice_dyn( 'V' ) 
    689677            ! 
    690678         ENDIF 
     
    701689      END DO                                              !  end loop over jter  ! 
    702690      !                                                   ! ==================== ! 
    703       ! 
    704       IF( ln_bdy .AND. .NOT.ll_bdy_substep ) THEN 
    705          CALL bdy_ice_dyn( 'U' ) 
    706          CALL bdy_ice_dyn( 'V' ) 
    707       ENDIF 
    708691      ! 
    709692      !------------------------------------------------------------------------------! 
     
    764747      DO jj = 1, jpj 
    765748         DO ji = 1, jpi 
    766             zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
     749            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
    767750         END DO 
    768751      END DO 
    769752 
     753      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 
     754      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 
     755         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 
     756         ! 
     757         CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., & 
     758            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 
     759         ! 
     760         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 
     761         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) 
     762         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 
     763         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 
     764         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 
     765         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 
     766      ENDIF 
     767        
    770768      ! --- divergence, shear and strength --- ! 
    771       IF( iom_use('icediv') )   CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:) )   ! divergence 
    772       IF( iom_use('iceshe') )   CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:) )   ! shear 
    773       IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength 
    774  
    775       ! --- charge ellipse --- ! 
    776       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN 
     769      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
     770      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
     771      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
     772 
     773      ! --- stress tensor --- ! 
     774      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
    777775         ! 
    778776         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
     
    780778         DO jj = 2, jpjm1 
    781779            DO ji = 2, jpim1 
    782                zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    783                   &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
    784                   &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 
     780               zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
     781                  &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
     782                  &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 
    785783 
    786784               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    787785 
    788                zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
     786               zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    789787 
    790788!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
     
    799797         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
    800798         ! 
    801          IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 ) 
    802          IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 ) 
    803          IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 ) 
    804          ! 
     799         CALL iom_put( 'isig1' , zsig1 ) 
     800         CALL iom_put( 'isig2' , zsig2 ) 
     801         CALL iom_put( 'isig3' , zsig3 ) 
     802         ! 
     803         ! Stress tensor invariants (normal and shear stress N/m) 
     804         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress 
     805         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
     806 
    805807         DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
    806808      ENDIF 
    807809       
    808810      ! --- SIMIP --- ! 
    809       IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. & 
    810          & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. & 
    811          & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. & 
    812          & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN 
    813  
    814          ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  & 
    815             &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  & 
    816             &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  & 
    817             &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) ) 
    818           
     811      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 
     812         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
     813         ! 
     814         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., & 
     815            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) 
     816 
     817         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x) 
     818         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y) 
     819         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x) 
     820         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y) 
     821         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x) 
     822         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y) 
     823      ENDIF 
     824 
     825      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. & 
     826         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN 
     827         ! 
     828         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & 
     829            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 
     830         ! 
    819831         DO jj = 2, jpjm1 
    820832            DO ji = 2, jpim1 
    821                rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
    822                 
    823                ! Stress tensor invariants (normal and shear stress N/m) 
    824                zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress 
    825                zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress 
    826                 
    827                ! Stress terms of the momentum equation (N/m2) 
    828                zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch     ! sea surface slope stress term 
    829                zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 
    830                 
    831                zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch     ! Coriolis stress term 
    832                zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch 
    833                 
    834                zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch     ! internal stress term 
    835                zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch 
    836                 
    837                zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress 
    838                zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 
    839                 
    840833               ! 2D ice mass, snow mass, area transport arrays (X, Y) 
    841                zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 
    842                zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
    843                 
     834               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 
     835               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 
     836 
    844837               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
    845838               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
    846                 
     839 
    847840               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
    848841               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
    849                 
     842 
    850843               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component 
    851844               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   '' 
    852                 
    853             END DO 
    854          END DO 
    855           
    856          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   & 
    857             &                zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   & 
    858             &                zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   &  
    859             &                zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    ) 
    860                    
    861          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_utau_oi  , 'U', -1., zdiag_vtau_oi  , 'V', -1.,   & 
    862             &                zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   & 
    863             &                zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   & 
    864             &                zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    ) 
    865           
    866          IF( iom_use('normstr' ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress 
    867          IF( iom_use('sheastr' ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress 
    868          IF( iom_use('dssh_dx' ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x) 
    869          IF( iom_use('dssh_dy' ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y) 
    870          IF( iom_use('corstrx' ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x) 
    871          IF( iom_use('corstry' ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y) 
    872          IF( iom_use('intstrx' ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x) 
    873          IF( iom_use('intstry' ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y) 
    874          IF( iom_use('utau_oi' ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x) 
    875          IF( iom_use('vtau_oi' ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y) 
    876          IF( iom_use('xmtrpice') )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s) 
    877          IF( iom_use('ymtrpice') )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport  
    878          IF( iom_use('xmtrpsnw') )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s) 
    879          IF( iom_use('ymtrpsnw') )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport 
    880          IF( iom_use('xatrp'   ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport 
    881          IF( iom_use('yatrp'   ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport 
    882  
    883          DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  & 
    884             &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  & 
    885             &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  & 
    886             &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     ) 
     845 
     846            END DO 
     847         END DO 
     848 
     849         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 
     850            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & 
     851            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. ) 
     852 
     853         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s) 
     854         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport  
     855         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s) 
     856         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport 
     857         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport 
     858         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport 
     859 
     860         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , & 
     861            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) 
    887862 
    888863      ENDIF 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icesbc.F90

    r11380 r11381  
    114114      INTEGER, INTENT(in) ::   ksbc   ! flux formulation (user defined, bulk or Pure Coupled) 
    115115      ! 
    116       INTEGER  ::   ji, jj, jl                                ! dummy loop index 
    117       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky 
    118       REAL(wp), DIMENSION(jpi,jpj)     ::   zalb              ! 2D workspace 
     116      INTEGER  ::   ji, jj, jl      ! dummy loop index 
     117      REAL(wp) ::   zmiss_val       ! missing value retrieved from xios  
     118      REAL(wp), DIMENSION(jpi,jpj,jpl)              ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky 
     119      REAL(wp), DIMENSION(:,:)        , ALLOCATABLE ::   zalb, zmsk00      ! 2D workspace 
    119120      !!-------------------------------------------------------------------- 
    120121      ! 
     
    126127         WRITE(numout,*)'~~~~~~~~~~~~~~~' 
    127128      ENDIF 
     129 
     130      ! get missing value from xml 
     131      CALL iom_miss_val( "icetemp", zmiss_val ) 
    128132 
    129133      ! --- cloud-sky and overcast-sky ice albedos --- ! 
     
    152156 
    153157      !--- output ice albedo and surface albedo ---! 
    154       IF( iom_use('icealb') ) THEN 
    155          WHERE( at_i_b <= epsi06 )   ;   zalb(:,:) = rn_alb_oce 
    156          ELSEWHERE                   ;   zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b 
     158      IF( iom_use('icealb') .OR. iom_use('albedo') ) THEN 
     159 
     160         ALLOCATE( zalb(jpi,jpj), zmsk00(jpi,jpj) ) 
     161 
     162         WHERE( at_i_b <= epsi06 ) 
     163            zmsk00(:,:) = 0._wp 
     164            zalb  (:,:) = rn_alb_oce 
     165         ELSEWHERE 
     166            zmsk00(:,:) = 1._wp             
     167            zalb  (:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b 
    157168         END WHERE 
    158          CALL iom_put( "icealb" , zalb(:,:) ) 
    159       ENDIF 
    160       IF( iom_use('albedo') ) THEN 
     169         ! ice albedo 
     170         CALL iom_put( 'icealb' , zalb * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) 
     171         ! ice+ocean albedo 
    161172         zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b ) 
    162          CALL iom_put( "albedo" , zalb(:,:) ) 
     173         CALL iom_put( 'albedo' , zalb ) 
     174 
     175         DEALLOCATE( zalb, zmsk00 ) 
     176 
    163177      ENDIF 
    164178      ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icestp.F90

    r11380 r11381  
    425425      wfx_err_sub(:,:) = 0._wp 
    426426      ! 
    427       afx_tot(:,:) = 0._wp   ; 
    428       ! 
    429427      diag_heat(:,:) = 0._wp ;   diag_sice(:,:) = 0._wp 
    430428      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/iceupdate.F90

    r11380 r11381  
    198198      ! --- salt fluxes [kg/m2/s] --- ! 
    199199      !                           ! sfxice =  sfxbog + sfxbom + sfxsum + sfxsni + sfxopw + sfxres + sfxdyn + sfxbri + sfxsub + sfxlam 
    200       IF( iom_use('sfxice'  ) )   CALL iom_put( "sfxice", sfx     * 1.e-03 )   ! salt flux from total ice growth/melt 
    201       IF( iom_use('sfxbog'  ) )   CALL iom_put( "sfxbog", sfx_bog * 1.e-03 )   ! salt flux from bottom growth 
    202       IF( iom_use('sfxbom'  ) )   CALL iom_put( "sfxbom", sfx_bom * 1.e-03 )   ! salt flux from bottom melting 
    203       IF( iom_use('sfxsum'  ) )   CALL iom_put( "sfxsum", sfx_sum * 1.e-03 )   ! salt flux from surface melting 
    204       IF( iom_use('sfxlam'  ) )   CALL iom_put( "sfxlam", sfx_lam * 1.e-03 )   ! salt flux from lateral melting 
    205       IF( iom_use('sfxsni'  ) )   CALL iom_put( "sfxsni", sfx_sni * 1.e-03 )   ! salt flux from snow ice formation 
    206       IF( iom_use('sfxopw'  ) )   CALL iom_put( "sfxopw", sfx_opw * 1.e-03 )   ! salt flux from open water formation 
    207       IF( iom_use('sfxdyn'  ) )   CALL iom_put( "sfxdyn", sfx_dyn * 1.e-03 )   ! salt flux from ridging rafting 
    208       IF( iom_use('sfxbri'  ) )   CALL iom_put( "sfxbri", sfx_bri * 1.e-03 )   ! salt flux from brines 
    209       IF( iom_use('sfxres'  ) )   CALL iom_put( "sfxres", sfx_res * 1.e-03 )   ! salt flux from undiagnosed processes 
    210       IF( iom_use('sfxsub'  ) )   CALL iom_put( "sfxsub", sfx_sub * 1.e-03 )   ! salt flux from sublimation 
     200      IF( iom_use('sfxice'  ) )   CALL iom_put( 'sfxice', sfx     * 1.e-03 )   ! salt flux from total ice growth/melt 
     201      IF( iom_use('sfxbog'  ) )   CALL iom_put( 'sfxbog', sfx_bog * 1.e-03 )   ! salt flux from bottom growth 
     202      IF( iom_use('sfxbom'  ) )   CALL iom_put( 'sfxbom', sfx_bom * 1.e-03 )   ! salt flux from bottom melting 
     203      IF( iom_use('sfxsum'  ) )   CALL iom_put( 'sfxsum', sfx_sum * 1.e-03 )   ! salt flux from surface melting 
     204      IF( iom_use('sfxlam'  ) )   CALL iom_put( 'sfxlam', sfx_lam * 1.e-03 )   ! salt flux from lateral melting 
     205      IF( iom_use('sfxsni'  ) )   CALL iom_put( 'sfxsni', sfx_sni * 1.e-03 )   ! salt flux from snow ice formation 
     206      IF( iom_use('sfxopw'  ) )   CALL iom_put( 'sfxopw', sfx_opw * 1.e-03 )   ! salt flux from open water formation 
     207      IF( iom_use('sfxdyn'  ) )   CALL iom_put( 'sfxdyn', sfx_dyn * 1.e-03 )   ! salt flux from ridging rafting 
     208      IF( iom_use('sfxbri'  ) )   CALL iom_put( 'sfxbri', sfx_bri * 1.e-03 )   ! salt flux from brines 
     209      IF( iom_use('sfxres'  ) )   CALL iom_put( 'sfxres', sfx_res * 1.e-03 )   ! salt flux from undiagnosed processes 
     210      IF( iom_use('sfxsub'  ) )   CALL iom_put( 'sfxsub', sfx_sub * 1.e-03 )   ! salt flux from sublimation 
    211211 
    212212      ! --- mass fluxes [kg/m2/s] --- ! 
    213       IF( iom_use('emp_oce' ) )   CALL iom_put( "emp_oce", emp_oce )   ! emp over ocean (taking into account the snow blown away from the ice) 
    214       IF( iom_use('emp_ice' ) )   CALL iom_put( "emp_ice", emp_ice )   ! emp over ice   (taking into account the snow blown away from the ice) 
     213      CALL iom_put( 'emp_oce', emp_oce )   ! emp over ocean (taking into account the snow blown away from the ice) 
     214      CALL iom_put( 'emp_ice', emp_ice )   ! emp over ice   (taking into account the snow blown away from the ice) 
    215215 
    216216      !                           ! vfxice = vfxbog + vfxbom + vfxsum + vfxsni + vfxopw + vfxdyn + vfxres + vfxlam + vfxpnd 
    217       IF( iom_use('vfxice'  ) )   CALL iom_put( "vfxice" , wfx_ice )   ! mass flux from total ice growth/melt 
    218       IF( iom_use('vfxbog'  ) )   CALL iom_put( "vfxbog" , wfx_bog )   ! mass flux from bottom growth 
    219       IF( iom_use('vfxbom'  ) )   CALL iom_put( "vfxbom" , wfx_bom )   ! mass flux from bottom melt  
    220       IF( iom_use('vfxsum'  ) )   CALL iom_put( "vfxsum" , wfx_sum )   ! mass flux from surface melt  
    221       IF( iom_use('vfxlam'  ) )   CALL iom_put( "vfxlam" , wfx_lam )   ! mass flux from lateral melt  
    222       IF( iom_use('vfxsni'  ) )   CALL iom_put( "vfxsni" , wfx_sni )   ! mass flux from snow-ice formation 
    223       IF( iom_use('vfxopw'  ) )   CALL iom_put( "vfxopw" , wfx_opw )   ! mass flux from growth in open water 
    224       IF( iom_use('vfxdyn'  ) )   CALL iom_put( "vfxdyn" , wfx_dyn )   ! mass flux from dynamics (ridging) 
    225       IF( iom_use('vfxres'  ) )   CALL iom_put( "vfxres" , wfx_res )   ! mass flux from undiagnosed processes  
    226       IF( iom_use('vfxpnd'  ) )   CALL iom_put( "vfxpnd" , wfx_pnd )   ! mass flux from melt ponds 
    227       IF( iom_use('vfxsub'  ) )   CALL iom_put( "vfxsub" , wfx_ice_sub )   ! mass flux from ice sublimation (ice-atm.) 
    228       IF( iom_use('vfxsub_err') ) CALL iom_put( "vfxsub_err", wfx_err_sub )   ! "excess" of sublimation sent to ocean       
    229  
    230       IF ( iom_use( "vfxthin" ) ) THEN   ! mass flux from ice growth in open water + thin ice (<20cm) => comparable to observations   
     217      CALL iom_put( 'vfxice'    , wfx_ice    )   ! mass flux from total ice growth/melt 
     218      CALL iom_put( 'vfxbog'    , wfx_bog    )   ! mass flux from bottom growth 
     219      CALL iom_put( 'vfxbom'    , wfx_bom    )   ! mass flux from bottom melt  
     220      CALL iom_put( 'vfxsum'    , wfx_sum    )   ! mass flux from surface melt  
     221      CALL iom_put( 'vfxlam'    , wfx_lam    )   ! mass flux from lateral melt  
     222      CALL iom_put( 'vfxsni'    , wfx_sni    )   ! mass flux from snow-ice formation 
     223      CALL iom_put( 'vfxopw'    , wfx_opw    )   ! mass flux from growth in open water 
     224      CALL iom_put( 'vfxdyn'    , wfx_dyn    )   ! mass flux from dynamics (ridging) 
     225      CALL iom_put( 'vfxres'    , wfx_res    )   ! mass flux from undiagnosed processes  
     226      CALL iom_put( 'vfxpnd'    , wfx_pnd    )   ! mass flux from melt ponds 
     227      CALL iom_put( 'vfxsub'    , wfx_ice_sub )   ! mass flux from ice sublimation (ice-atm.) 
     228      CALL iom_put( 'vfxsub_err', wfx_err_sub )   ! "excess" of sublimation sent to ocean       
     229 
     230      IF ( iom_use( 'vfxthin' ) ) THEN   ! mass flux from ice growth in open water + thin ice (<20cm) => comparable to observations   
    231231         WHERE( hm_i(:,:) < 0.2 .AND. hm_i(:,:) > 0. ) ; z2d = wfx_bog 
    232232         ELSEWHERE                                     ; z2d = 0._wp 
    233233         END WHERE 
    234          CALL iom_put( "vfxthin", wfx_opw + z2d ) 
    235       ENDIF 
    236  
    237       !                              ! vfxsnw = vfxsnw_sni + vfxsnw_dyn + vfxsnw_sum 
    238       IF( iom_use('vfxsnw'     ) )   CALL iom_put( "vfxsnw"     , wfx_snw     )   ! mass flux from total snow growth/melt 
    239       IF( iom_use('vfxsnw_sum' ) )   CALL iom_put( "vfxsnw_sum" , wfx_snw_sum )   ! mass flux from snow melt at the surface 
    240       IF( iom_use('vfxsnw_sni' ) )   CALL iom_put( "vfxsnw_sni" , wfx_snw_sni )   ! mass flux from snow melt during snow-ice formation  
    241       IF( iom_use('vfxsnw_dyn' ) )   CALL iom_put( "vfxsnw_dyn" , wfx_snw_dyn )   ! mass flux from dynamics (ridging)  
    242       IF( iom_use('vfxsnw_sub' ) )   CALL iom_put( "vfxsnw_sub" , wfx_snw_sub )   ! mass flux from snow sublimation (ice-atm.)  
    243       IF( iom_use('vfxsnw_pre' ) )   CALL iom_put( "vfxsnw_pre" , wfx_spr     )   ! snow precip 
     234         CALL iom_put( 'vfxthin', wfx_opw + z2d ) 
     235      ENDIF 
     236 
     237      !                            ! vfxsnw = vfxsnw_sni + vfxsnw_dyn + vfxsnw_sum 
     238      CALL iom_put( 'vfxsnw'     , wfx_snw     )   ! mass flux from total snow growth/melt 
     239      CALL iom_put( 'vfxsnw_sum' , wfx_snw_sum )   ! mass flux from snow melt at the surface 
     240      CALL iom_put( 'vfxsnw_sni' , wfx_snw_sni )   ! mass flux from snow melt during snow-ice formation  
     241      CALL iom_put( 'vfxsnw_dyn' , wfx_snw_dyn )   ! mass flux from dynamics (ridging)  
     242      CALL iom_put( 'vfxsnw_sub' , wfx_snw_sub )   ! mass flux from snow sublimation (ice-atm.)  
     243      CALL iom_put( 'vfxsnw_pre' , wfx_spr     )   ! snow precip 
    244244 
    245245      ! --- heat fluxes [W/m2] --- ! 
    246246      !                              ! qt_atm_oi - qt_oce_ai = hfxdhc - ( dihctrp + dshctrp ) 
    247       IF( iom_use('qsr_oce'    ) )   CALL iom_put( "qsr_oce"    , qsr_oce * ( 1._wp - at_i_b )                               )   !     solar flux at ocean surface 
    248       IF( iom_use('qns_oce'    ) )   CALL iom_put( "qns_oce"    , qns_oce * ( 1._wp - at_i_b ) + qemp_oce                    )   ! non-solar flux at ocean surface 
    249       IF( iom_use('qsr_ice'    ) )   CALL iom_put( "qsr_ice"    , SUM( qsr_ice * a_i_b, dim=3 )                              )   !     solar flux at ice surface 
    250       IF( iom_use('qns_ice'    ) )   CALL iom_put( "qns_ice"    , SUM( qns_ice * a_i_b, dim=3 ) + qemp_ice                   )   ! non-solar flux at ice surface 
    251       IF( iom_use('qtr_ice_bot') )   CALL iom_put( "qtr_ice_bot", SUM( qtr_ice_bot * a_i_b, dim=3 )                          )   !     solar flux transmitted thru ice 
    252       IF( iom_use('qtr_ice_top') )   CALL iom_put( "qtr_ice_top", SUM( qtr_ice_top * a_i_b, dim=3 )                          )   !     solar flux transmitted thru ice surface 
    253       IF( iom_use('qt_oce'     ) )   CALL iom_put( "qt_oce"     ,      ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce ) 
    254       IF( iom_use('qt_ice'     ) )   CALL iom_put( "qt_ice"     , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 )     + qemp_ice ) 
    255       IF( iom_use('qt_oce_ai'  ) )   CALL iom_put( "qt_oce_ai"  , qt_oce_ai * tmask(:,:,1)                                   )   ! total heat flux at the ocean   surface: interface oce-(ice+atm)  
    256       IF( iom_use('qt_atm_oi'  ) )   CALL iom_put( "qt_atm_oi"  , qt_atm_oi * tmask(:,:,1)                                   )   ! total heat flux at the oce-ice surface: interface atm-(ice+oce)  
    257       IF( iom_use('qemp_oce'   ) )   CALL iom_put( "qemp_oce"   , qemp_oce                                                   )   ! Downward Heat Flux from E-P over ocean 
    258       IF( iom_use('qemp_ice'   ) )   CALL iom_put( "qemp_ice"   , qemp_ice                                                   )   ! Downward Heat Flux from E-P over ice 
     247      IF( iom_use('qsr_oce'    ) )   CALL iom_put( 'qsr_oce'    , qsr_oce * ( 1._wp - at_i_b )                               )   !     solar flux at ocean surface 
     248      IF( iom_use('qns_oce'    ) )   CALL iom_put( 'qns_oce'    , qns_oce * ( 1._wp - at_i_b ) + qemp_oce                    )   ! non-solar flux at ocean surface 
     249      IF( iom_use('qsr_ice'    ) )   CALL iom_put( 'qsr_ice'    , SUM( qsr_ice * a_i_b, dim=3 )                              )   !     solar flux at ice surface 
     250      IF( iom_use('qns_ice'    ) )   CALL iom_put( 'qns_ice'    , SUM( qns_ice * a_i_b, dim=3 ) + qemp_ice                   )   ! non-solar flux at ice surface 
     251      IF( iom_use('qtr_ice_bot') )   CALL iom_put( 'qtr_ice_bot', SUM( qtr_ice_bot * a_i_b, dim=3 )                          )   !     solar flux transmitted thru ice 
     252      IF( iom_use('qtr_ice_top') )   CALL iom_put( 'qtr_ice_top', SUM( qtr_ice_top * a_i_b, dim=3 )                          )   !     solar flux transmitted thru ice surface 
     253      IF( iom_use('qt_oce'     ) )   CALL iom_put( 'qt_oce'     ,      ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce ) 
     254      IF( iom_use('qt_ice'     ) )   CALL iom_put( 'qt_ice'     , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 )     + qemp_ice ) 
     255      IF( iom_use('qt_oce_ai'  ) )   CALL iom_put( 'qt_oce_ai'  , qt_oce_ai * tmask(:,:,1)                                   )   ! total heat flux at the ocean   surface: interface oce-(ice+atm)  
     256      IF( iom_use('qt_atm_oi'  ) )   CALL iom_put( 'qt_atm_oi'  , qt_atm_oi * tmask(:,:,1)                                   )   ! total heat flux at the oce-ice surface: interface atm-(ice+oce)  
     257      IF( iom_use('qemp_oce'   ) )   CALL iom_put( 'qemp_oce'   , qemp_oce                                                   )   ! Downward Heat Flux from E-P over ocean 
     258      IF( iom_use('qemp_ice'   ) )   CALL iom_put( 'qemp_ice'   , qemp_ice                                                   )   ! Downward Heat Flux from E-P over ice 
    259259 
    260260      ! heat fluxes from ice transformations 
    261       !                              ! hfxdhc = hfxbog + hfxbom + hfxsum + hfxopw + hfxdif + hfxsnw - ( hfxthd + hfxdyn + hfxres + hfxsub + hfxspr ) 
    262       IF( iom_use('hfxbog'     ) )   CALL iom_put ("hfxbog"     , hfx_bog             )   ! heat flux used for ice bottom growth  
    263       IF( iom_use('hfxbom'     ) )   CALL iom_put ("hfxbom"     , hfx_bom             )   ! heat flux used for ice bottom melt 
    264       IF( iom_use('hfxsum'     ) )   CALL iom_put ("hfxsum"     , hfx_sum             )   ! heat flux used for ice surface melt 
    265       IF( iom_use('hfxopw'     ) )   CALL iom_put ("hfxopw"     , hfx_opw             )   ! heat flux used for ice formation in open water 
    266       IF( iom_use('hfxdif'     ) )   CALL iom_put ("hfxdif"     , hfx_dif             )   ! heat flux used for ice temperature change 
    267       IF( iom_use('hfxsnw'     ) )   CALL iom_put ("hfxsnw"     , hfx_snw             )   ! heat flux used for snow melt  
    268       IF( iom_use('hfxerr'     ) )   CALL iom_put ("hfxerr"     , hfx_err_dif        )   ! heat flux error after heat diffusion (included in qt_oce_ai) 
     261      !                            ! hfxdhc = hfxbog + hfxbom + hfxsum + hfxopw + hfxdif + hfxsnw - ( hfxthd + hfxdyn + hfxres + hfxsub + hfxspr ) 
     262      CALL iom_put ('hfxbog'     , hfx_bog     )   ! heat flux used for ice bottom growth  
     263      CALL iom_put ('hfxbom'     , hfx_bom     )   ! heat flux used for ice bottom melt 
     264      CALL iom_put ('hfxsum'     , hfx_sum     )   ! heat flux used for ice surface melt 
     265      CALL iom_put ('hfxopw'     , hfx_opw     )   ! heat flux used for ice formation in open water 
     266      CALL iom_put ('hfxdif'     , hfx_dif     )   ! heat flux used for ice temperature change 
     267      CALL iom_put ('hfxsnw'     , hfx_snw     )   ! heat flux used for snow melt  
     268      CALL iom_put ('hfxerr'     , hfx_err_dif )   ! heat flux error after heat diffusion (included in qt_oce_ai) 
    269269 
    270270      ! heat fluxes associated with mass exchange (freeze/melt/precip...) 
    271       IF( iom_use('hfxthd'     ) )   CALL iom_put ("hfxthd"     , hfx_thd             )   !   
    272       IF( iom_use('hfxdyn'     ) )   CALL iom_put ("hfxdyn"     , hfx_dyn             )   !   
    273       IF( iom_use('hfxres'     ) )   CALL iom_put ("hfxres"     , hfx_res             )   !   
    274       IF( iom_use('hfxsub'     ) )   CALL iom_put ("hfxsub"     , hfx_sub             )   !   
    275       IF( iom_use('hfxspr'     ) )   CALL iom_put ("hfxspr"     , hfx_spr             )   ! Heat flux from snow precip heat content  
     271      CALL iom_put ('hfxthd'     , hfx_thd     )   !   
     272      CALL iom_put ('hfxdyn'     , hfx_dyn     )   !   
     273      CALL iom_put ('hfxres'     , hfx_res     )   !   
     274      CALL iom_put ('hfxsub'     , hfx_sub     )   !   
     275      CALL iom_put ('hfxspr'     , hfx_spr     )   ! Heat flux from snow precip heat content  
    276276 
    277277      ! other heat fluxes 
    278       IF( iom_use('hfxsensib'  ) )   CALL iom_put( "hfxsensib"  ,     -qsb_ice_bot * at_i_b         )   ! Sensible oceanic heat flux 
    279       IF( iom_use('hfxcndbot'  ) )   CALL iom_put( "hfxcndbot"  , SUM( qcn_ice_bot * a_i_b, dim=3 ) )   ! Bottom conduction flux 
    280       IF( iom_use('hfxcndtop'  ) )   CALL iom_put( "hfxcndtop"  , SUM( qcn_ice_top * a_i_b, dim=3 ) )   ! Surface conduction flux 
     278      IF( iom_use('hfxsensib'  ) )   CALL iom_put( 'hfxsensib'  ,     -qsb_ice_bot * at_i_b         )   ! Sensible oceanic heat flux 
     279      IF( iom_use('hfxcndbot'  ) )   CALL iom_put( 'hfxcndbot'  , SUM( qcn_ice_bot * a_i_b, dim=3 ) )   ! Bottom conduction flux 
     280      IF( iom_use('hfxcndtop'  ) )   CALL iom_put( 'hfxcndtop'  , SUM( qcn_ice_top * a_i_b, dim=3 ) )   ! Surface conduction flux 
    281281 
    282282      ! diags 
    283       IF( iom_use('hfxdhc'     ) )   CALL iom_put ("hfxdhc"     , diag_heat           )   ! Heat content variation in snow and ice  
     283      CALL iom_put ('hfxdhc'     , diag_heat   )   ! Heat content variation in snow and ice  
    284284      ! 
    285285      ! controls 
     
    413413      !! ** Method  :   use of IOM library 
    414414      !!---------------------------------------------------------------------- 
    415       CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     415      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! 'READ'/'WRITE' flag 
    416416      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step 
    417417      ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icevar.F90

    r11380 r11381  
    3232   !!                        - vt_s(jpi,jpj) 
    3333   !!                        - at_i(jpi,jpj) 
     34   !!                        - st_i(jpi,jpj) 
    3435   !!                        - et_s(jpi,jpj)  total snow heat content 
    3536   !!                        - et_i(jpi,jpj)  total ice thermal content  
     
    104105      ! 
    105106      !                                      ! integrated values 
    106       vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 ) 
    107       vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 ) 
    108       at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 ) 
    109       et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 
    110       et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 
     107      vt_i(:,:) =       SUM( v_i (:,:,:)           , dim=3 ) 
     108      vt_s(:,:) =       SUM( v_s (:,:,:)           , dim=3 ) 
     109      st_i(:,:) =       SUM( sv_i(:,:,:)           , dim=3 ) 
     110      at_i(:,:) =       SUM( a_i (:,:,:)           , dim=3 ) 
     111      et_s(:,:)  = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 ) 
     112      et_i(:,:)  = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 ) 
    111113      ! 
    112114      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds 
     
    138140         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    139141         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:) 
    140          sm_i (:,:) = SUM( sv_i(:,:,:)              , dim=3 ) * z1_vt_i(:,:) 
     142         sm_i (:,:) =      st_i(:,:)                          * z1_vt_i(:,:) 
    141143         ! 
    142144         tm_i(:,:) = 0._wp 
     
    263265      ! 
    264266      ! integrated values  
    265       vt_i (:,:) = SUM( v_i, dim=3 ) 
    266       vt_s (:,:) = SUM( v_s, dim=3 ) 
    267       at_i (:,:) = SUM( a_i, dim=3 ) 
     267      vt_i (:,:) = SUM( v_i , dim=3 ) 
     268      vt_s (:,:) = SUM( v_s , dim=3 ) 
     269      at_i (:,:) = SUM( a_i , dim=3 ) 
    268270      ! 
    269271   END SUBROUTINE ice_var_glo2eqv 
     
    533535 
    534536      ! to be sure that at_i is the sum of a_i(jl) 
    535       at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) 
    536       vt_i (:,:) = SUM( v_i(:,:,:), dim=3 ) 
     537      at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) 
     538      vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) 
     539!!clem add? 
     540!      vt_s (:,:) = SUM( v_s (:,:,:), dim=3 ) 
     541!      st_i (:,:) = SUM( sv_i(:,:,:), dim=3 ) 
     542!      et_s(:,:)  = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 ) 
     543!      et_i(:,:)  = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 ) 
     544!!clem 
    537545 
    538546      ! open water = 1 if at_i=0 
     
    10851093         !                              ! ---------------------- ! 
    10861094         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:) ) 
    1087 !!$         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:), & 
    1088 !!$            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) ) 
     1095!!         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:), & 
     1096!!            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) ) 
    10891097         !                              ! ---------------------- ! 
    10901098      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         ! 
    10911099         !                              ! ---------------------- ! 
    1092          CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1) )          
    1093 !!$         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1), & 
    1094 !!$            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) )          
     1100         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1) ) 
     1101!!         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1), & 
     1102!!            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) )          
    10951103         !                              ! ----------------------- ! 
    10961104      ELSE                              ! input cat /= output cat ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icewri.F90

    r11380 r11381  
    5050      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices 
    5151      REAL(wp) ::   z2da, z2db, zrho1, zrho2 
    52       REAL(wp), DIMENSION(jpi,jpj)     ::   z2d, zfast !  2D workspace 
     52      REAL(wp) ::   zmiss_val       ! missing value retrieved from xios  
     53      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d, zfast                     ! 2D workspace 
    5354      REAL(wp), DIMENSION(jpi,jpj)     ::   zmsk00, zmsk05, zmsk15, zmsksn ! O%, 5% and 15% concentration mask and snow mask 
    5455      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zmsk00l, zmsksnl               ! cat masks 
     
    5859      REAL(wp) ::   zdiag_area_sh, zdiag_extt_sh, zdiag_volu_sh  
    5960      !!------------------------------------------------------------------- 
    60  
     61      ! 
    6162      IF( ln_timing )   CALL timing_start('icewri') 
     63 
     64      ! get missing value from xml 
     65      CALL iom_miss_val( 'icetemp', zmiss_val ) 
    6266 
    6367      ! brine volume 
     
    8589      ! Standard outputs 
    8690      !----------------- 
    87       zrho1 = ( rau0 - rhoi ) * r1_rau0; zrho2 = rhos * r1_rau0 
     91      zrho1 = ( rau0 - rhoi ) * r1_rau0 ; zrho2 = rhos * r1_rau0 
    8892      ! masks 
    89       IF( iom_use('icemask'  ) )   CALL iom_put( "icemask"  , zmsk00              )   ! ice mask 0% 
    90       IF( iom_use('icemask05') )   CALL iom_put( "icemask05", zmsk05              )   ! ice mask 5% 
    91       IF( iom_use('icemask15') )   CALL iom_put( "icemask15", zmsk15              )   ! ice mask 15% 
     93      CALL iom_put( 'icemask'  , zmsk00 )   ! ice mask 0% 
     94      CALL iom_put( 'icemask05', zmsk05 )   ! ice mask 5% 
     95      CALL iom_put( 'icemask15', zmsk15 )   ! ice mask 15% 
     96      CALL iom_put( 'icepres'  , zmsk00 )   ! Ice presence (1 or 0)  
    9297      ! 
    9398      ! general fields 
    94       IF( iom_use('icemass'  ) )   CALL iom_put( "icemass", rhoi * vt_i * zmsk00  )   ! Ice mass per cell area  
    95       IF( iom_use('snwmass'  ) )   CALL iom_put( "snwmass", rhos * vt_s * zmsksn  )   ! Snow mass per cell area 
    96       IF( iom_use('icepres'  ) )   CALL iom_put( "icepres", zmsk00                )   ! Ice presence (1 or 0)  
    97       IF( iom_use('iceconc'  ) )   CALL iom_put( "iceconc", at_i  * zmsk00        )   ! ice concentration 
    98       IF( iom_use('icevolu'  ) )   CALL iom_put( "icevolu", vt_i  * zmsk00        )   ! ice volume = mean ice thickness over the cell 
    99       IF( iom_use('icethic'  ) )   CALL iom_put( "icethic", hm_i  * zmsk00        )   ! ice thickness 
    100       IF( iom_use('snwthic'  ) )   CALL iom_put( "snwthic", hm_s  * zmsk00        )   ! snw thickness 
    101       IF( iom_use('icebrv'   ) )   CALL iom_put( "icebrv" , bvm_i * zmsk00 * 100. )   ! brine volume 
    102       IF( iom_use('iceage'   ) )   CALL iom_put( "iceage" , om_i  * zmsk15 / rday )   ! ice age 
    103       IF( iom_use('icehnew'  ) )   CALL iom_put( "icehnew", ht_i_new              )   ! new ice thickness formed in the leads 
    104       IF( iom_use('snwvolu'  ) )   CALL iom_put( "snwvolu", vt_s  * zmsksn        )   ! snow volume 
    105       IF( iom_use('icefrb') ) THEN 
     99      IF( iom_use('icemass' ) )   CALL iom_put( 'icemass', vt_i * rhoi * zmsk00 )                                           ! Ice mass per cell area  
     100      IF( iom_use('snwmass' ) )   CALL iom_put( 'snwmass', vt_s * rhos * zmsksn )                                           ! Snow mass per cell area 
     101      IF( iom_use('iceconc' ) )   CALL iom_put( 'iceconc', at_i        * zmsk00 )                                           ! ice concentration 
     102      IF( iom_use('icevolu' ) )   CALL iom_put( 'icevolu', vt_i        * zmsk00 )                                           ! ice volume = mean ice thickness over the cell 
     103      IF( iom_use('icethic' ) )   CALL iom_put( 'icethic', hm_i        * zmsk00 )                                           ! ice thickness 
     104      IF( iom_use('snwthic' ) )   CALL iom_put( 'snwthic', hm_s        * zmsk00 )                                           ! snw thickness 
     105      IF( iom_use('icebrv'  ) )   CALL iom_put( 'icebrv' , bvm_i* 100. * zmsk00 )                                           ! brine volume 
     106      IF( iom_use('iceage'  ) )   CALL iom_put( 'iceage' , om_i / rday * zmsk15 + zmiss_val * ( 1._wp - zmsk15 ) )          ! ice age 
     107      IF( iom_use('icehnew' ) )   CALL iom_put( 'icehnew', ht_i_new             )                                           ! new ice thickness formed in the leads 
     108      IF( iom_use('snwvolu' ) )   CALL iom_put( 'snwvolu', vt_s        * zmsksn )                                           ! snow volume 
     109      IF( iom_use('icefrb'  ) ) THEN                                                                                        ! Ice freeboard 
    106110         z2d(:,:) = ( zrho1 * hm_i(:,:) - zrho2 * hm_s(:,:) )                                          
    107111         WHERE( z2d < 0._wp )   z2d = 0._wp 
    108                                    CALL iom_put( "icefrb" , z2d * zmsk00          )   ! Ice freeboard 
     112                                  CALL iom_put( 'icefrb' , z2d * zmsk00         ) 
    109113      ENDIF 
    110       ! 
    111114      ! melt ponds 
    112       IF( iom_use('iceapnd'  ) )   CALL iom_put( "iceapnd", at_ip  * zmsk00       )   ! melt pond total fraction 
    113       IF( iom_use('icevpnd'  ) )   CALL iom_put( "icevpnd", vt_ip  * zmsk00       )   ! melt pond total volume per unit area 
    114       ! 
     115      IF( iom_use('iceapnd' ) )   CALL iom_put( 'iceapnd', at_ip  * zmsk00      )                                           ! melt pond total fraction 
     116      IF( iom_use('icevpnd' ) )   CALL iom_put( 'icevpnd', vt_ip  * zmsk00      )                                           ! melt pond total volume per unit area 
    115117      ! salt 
    116       IF( iom_use('icesalt'  ) )   CALL iom_put( "icesalt", sm_i  * zmsk00        )   ! mean ice salinity 
    117       IF( iom_use('icesalm'  ) )   CALL iom_put( "icesalm", SUM( sv_i, DIM = 3 ) * rhoi * 1.0e-3 * zmsk00 )   ! Mass of salt in sea ice per cell area 
    118  
     118      IF( iom_use('icesalt' ) )   CALL iom_put( 'icesalt', sm_i                 * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! mean ice salinity 
     119      IF( iom_use('icesalm' ) )   CALL iom_put( 'icesalm', st_i * rhoi * 1.0e-3 * zmsk00 )                                  ! Mass of salt in sea ice per cell area 
    119120      ! heat 
    120       IF( iom_use('icetemp'  ) )   CALL iom_put( "icetemp", ( tm_i  - rt0 ) * zmsk00 )   ! ice mean temperature 
    121       IF( iom_use('snwtemp'  ) )   CALL iom_put( "snwtemp", ( tm_s  - rt0 ) * zmsksn )   ! snw mean temperature 
    122       IF( iom_use('icettop'  ) )   CALL iom_put( "icettop", ( tm_su - rt0 ) * zmsk00 )   ! temperature at the ice surface 
    123       IF( iom_use('icetbot'  ) )   CALL iom_put( "icetbot", ( t_bo  - rt0 ) * zmsk00 )   ! temperature at the ice bottom 
    124       IF( iom_use('icetsni'  ) )   CALL iom_put( "icetsni", ( tm_si - rt0 ) * zmsk00 )   ! temperature at the snow-ice interface 
    125       IF( iom_use('icehc'    ) )   CALL iom_put( "icehc"  ,  -et_i          * zmsk00 )   ! ice heat content 
    126       IF( iom_use('snwhc'    ) )   CALL iom_put( "snwhc"  ,  -et_s          * zmsksn )   ! snow heat content 
    127  
     121      IF( iom_use('icetemp' ) )   CALL iom_put( 'icetemp', ( tm_i  - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! ice mean temperature 
     122      IF( iom_use('snwtemp' ) )   CALL iom_put( 'snwtemp', ( tm_s  - rt0 ) * zmsksn + zmiss_val * ( 1._wp - zmsksn ) )      ! snw mean temperature 
     123      IF( iom_use('icettop' ) )   CALL iom_put( 'icettop', ( tm_su - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the ice surface 
     124      IF( iom_use('icetbot' ) )   CALL iom_put( 'icetbot', ( t_bo  - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the ice bottom 
     125      IF( iom_use('icetsni' ) )   CALL iom_put( 'icetsni', ( tm_si - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the snow-ice interface 
     126      IF( iom_use('icehc'   ) )   CALL iom_put( 'icehc'  ,  -et_i          * zmsk00 )                                       ! ice heat content 
     127      IF( iom_use('snwhc'   ) )   CALL iom_put( 'snwhc'  ,  -et_s          * zmsksn )                                       ! snow heat content 
    128128      ! momentum 
    129       IF( iom_use('uice'     ) )   CALL iom_put( "uice"   , u_ice                 )   ! ice velocity u component 
    130       IF( iom_use('vice'     ) )   CALL iom_put( "vice"   , v_ice                 )   ! ice velocity v component 
    131       IF( iom_use('utau_ai'  ) )   CALL iom_put( "utau_ai", utau_ice * zmsk00     )   ! Wind stress term in force balance (x) 
    132       IF( iom_use('vtau_ai'  ) )   CALL iom_put( "vtau_ai", vtau_ice * zmsk00     )   ! Wind stress term in force balance (y) 
    133  
    134       IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN  
    135         ! module of ice velocity 
     129      IF( iom_use('uice'    ) )   CALL iom_put( 'uice'   , u_ice    )                                                       ! ice velocity u 
     130      IF( iom_use('vice'    ) )   CALL iom_put( 'vice'   , v_ice    )                                                       ! ice velocity v 
     131      ! 
     132      IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN                                                              ! module of ice velocity 
    136133         DO jj = 2 , jpjm1 
    137134            DO ji = 2 , jpim1 
    138                z2da  = ( u_ice(ji,jj) + u_ice(ji-1,jj) ) 
    139                z2db  = ( v_ice(ji,jj) + v_ice(ji,jj-1) ) 
     135               z2da  = u_ice(ji,jj) + u_ice(ji-1,jj) 
     136               z2db  = v_ice(ji,jj) + v_ice(ji,jj-1) 
    140137               z2d(ji,jj) = 0.5_wp * SQRT( z2da * z2da + z2db * z2db ) 
    141138           END DO 
    142139         END DO 
    143140         CALL lbc_lnk( 'icewri', z2d, 'T', 1. ) 
    144          IF( iom_use('icevel') )   CALL iom_put( "icevel" , z2d ) 
    145  
    146         ! record presence of fast ice 
    147          WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp 
     141         CALL iom_put( 'icevel', z2d ) 
     142 
     143         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp                                      ! record presence of fast ice 
    148144         ELSEWHERE                                                ; zfast(:,:) = 0._wp 
    149145         END WHERE 
    150          IF( iom_use('fasticepres') )   CALL iom_put( "fasticepres" , zfast ) 
     146         CALL iom_put( 'fasticepres', zfast ) 
    151147      ENDIF 
    152148 
    153149      ! --- category-dependent fields --- ! 
    154       IF( iom_use('icemask_cat' ) )   CALL iom_put( "icemask_cat" , zmsk00l                                                    )   ! ice mask 0% 
    155       IF( iom_use('iceconc_cat' ) )   CALL iom_put( "iceconc_cat" , a_i * zmsk00l                                              )   ! area for categories 
    156       IF( iom_use('icethic_cat' ) )   CALL iom_put( "icethic_cat" , h_i * zmsk00l                                              )   ! thickness for categories 
    157       IF( iom_use('snwthic_cat' ) )   CALL iom_put( "snwthic_cat" , h_s * zmsksnl                                              )   ! snow depth for categories 
    158       IF( iom_use('icesalt_cat' ) )   CALL iom_put( "icesalt_cat" , s_i * zmsk00l                                              )   ! salinity for categories 
    159       IF( iom_use('iceage_cat'  ) )   CALL iom_put( "iceage_cat"  , o_i * zmsk00l / rday                                       )   ! ice age 
    160       IF( iom_use('icetemp_cat' ) )   CALL iom_put( "icetemp_cat" , ( SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i - rt0 ) * zmsk00l )   ! ice temperature 
    161       IF( iom_use('snwtemp_cat' ) )   CALL iom_put( "snwtemp_cat" , ( SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s - rt0 ) * zmsksnl )   ! snow temperature 
    162       IF( iom_use('icettop_cat' ) )   CALL iom_put( "icettop_cat" , ( t_su - rt0 ) * zmsk00l                                   )   ! surface temperature 
    163       IF( iom_use('icebrv_cat'  ) )   CALL iom_put( "icebrv_cat"  ,   bv_i * 100.  * zmsk00l                                   )   ! brine volume 
    164       IF( iom_use('iceapnd_cat' ) )   CALL iom_put( "iceapnd_cat" ,   a_ip         * zmsk00l                                   )   ! melt pond frac for categories 
    165       IF( iom_use('icehpnd_cat' ) )   CALL iom_put( "icehpnd_cat" ,   h_ip         * zmsk00l                                   )   ! melt pond frac for categories 
    166       IF( iom_use('iceafpnd_cat') )   CALL iom_put( "iceafpnd_cat",   a_ip_frac    * zmsk00l                                   )   ! melt pond frac for categories 
     150      IF( iom_use('icemask_cat' ) )   CALL iom_put( 'icemask_cat' ,                  zmsk00l                                   ) ! ice mask 0% 
     151      IF( iom_use('iceconc_cat' ) )   CALL iom_put( 'iceconc_cat' , a_i            * zmsk00l                                   ) ! area for categories 
     152      IF( iom_use('icethic_cat' ) )   CALL iom_put( 'icethic_cat' , h_i            * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! thickness for categories 
     153      IF( iom_use('snwthic_cat' ) )   CALL iom_put( 'snwthic_cat' , h_s            * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow depth for categories 
     154      IF( iom_use('icesalt_cat' ) )   CALL iom_put( 'icesalt_cat' , s_i            * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! salinity for categories 
     155      IF( iom_use('iceage_cat'  ) )   CALL iom_put( 'iceage_cat'  , o_i / rday     * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice age 
     156      IF( iom_use('icetemp_cat' ) )   CALL iom_put( 'icetemp_cat' , ( SUM( t_i, dim=3 ) * r1_nlay_i - rt0 ) & 
     157         &                                                                         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice temperature 
     158      IF( iom_use('snwtemp_cat' ) )   CALL iom_put( 'snwtemp_cat' , ( SUM( t_s, dim=3 ) * r1_nlay_s - rt0 ) & 
     159         &                                                                         * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow temperature 
     160      IF( iom_use('icettop_cat' ) )   CALL iom_put( 'icettop_cat' , ( t_su - rt0 ) * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! surface temperature 
     161      IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume 
     162      IF( iom_use('iceapnd_cat' ) )   CALL iom_put( 'iceapnd_cat' ,   a_ip         * zmsk00l                                   ) ! melt pond frac for categories 
     163      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond frac for categories 
     164      IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories 
    167165 
    168166      !------------------ 
     
    170168      !------------------ 
    171169      ! trends 
    172       IF( iom_use('dmithd') )   CALL iom_put( "dmithd", - wfx_bog - wfx_bom - wfx_sum - wfx_sni - wfx_opw - wfx_lam - wfx_res ) ! Sea-ice mass change from thermodynamics 
    173       IF( iom_use('dmidyn') )   CALL iom_put( "dmidyn", - wfx_dyn + rhoi * diag_trp_vi      )  ! Sea-ice mass change from dynamics(kg/m2/s) 
    174       IF( iom_use('dmiopw') )   CALL iom_put( "dmiopw", - wfx_opw                           )  ! Sea-ice mass change through growth in open water 
    175       IF( iom_use('dmibog') )   CALL iom_put( "dmibog", - wfx_bog                           )  ! Sea-ice mass change through basal growth 
    176       IF( iom_use('dmisni') )   CALL iom_put( "dmisni", - wfx_sni                           )  ! Sea-ice mass change through snow-to-ice conversion 
    177       IF( iom_use('dmisum') )   CALL iom_put( "dmisum", - wfx_sum                           )  ! Sea-ice mass change through surface melting 
    178       IF( iom_use('dmibom') )   CALL iom_put( "dmibom", - wfx_bom                           )  ! Sea-ice mass change through bottom melting 
    179       IF( iom_use('dmtsub') )   CALL iom_put( "dmtsub", - wfx_sub                           )  ! Sea-ice mass change through evaporation and sublimation 
    180       IF( iom_use('dmssub') )   CALL iom_put( "dmssub", - wfx_snw_sub                       )  ! Snow mass change through sublimation 
    181       IF( iom_use('dmisub') )   CALL iom_put( "dmisub", - wfx_ice_sub                       )  ! Sea-ice mass change through sublimation 
    182       IF( iom_use('dmsspr') )   CALL iom_put( "dmsspr", - wfx_spr                           )  ! Snow mass change through snow fall 
    183       IF( iom_use('dmsssi') )   CALL iom_put( "dmsssi",   wfx_sni*rhos*r1_rhoi              )  ! Snow mass change through snow-to-ice conversion 
    184       IF( iom_use('dmsmel') )   CALL iom_put( "dmsmel", - wfx_snw_sum                       )  ! Snow mass change through melt 
    185       IF( iom_use('dmsdyn') )   CALL iom_put( "dmsdyn", - wfx_snw_dyn + rhos * diag_trp_vs  )  ! Snow mass change through dynamics(kg/m2/s) 
    186  
     170      IF( iom_use('dmithd') )   CALL iom_put( 'dmithd', - wfx_bog - wfx_bom - wfx_sum - wfx_sni - wfx_opw - wfx_lam - wfx_res ) ! Sea-ice mass change from thermodynamics 
     171      IF( iom_use('dmidyn') )   CALL iom_put( 'dmidyn', - wfx_dyn + rhoi * diag_trp_vi                                        ) ! Sea-ice mass change from dynamics(kg/m2/s) 
     172      IF( iom_use('dmiopw') )   CALL iom_put( 'dmiopw', - wfx_opw                                                             ) ! Sea-ice mass change through growth in open water 
     173      IF( iom_use('dmibog') )   CALL iom_put( 'dmibog', - wfx_bog                                                             ) ! Sea-ice mass change through basal growth 
     174      IF( iom_use('dmisni') )   CALL iom_put( 'dmisni', - wfx_sni                                                             ) ! Sea-ice mass change through snow-to-ice conversion 
     175      IF( iom_use('dmisum') )   CALL iom_put( 'dmisum', - wfx_sum                                                             ) ! Sea-ice mass change through surface melting 
     176      IF( iom_use('dmibom') )   CALL iom_put( 'dmibom', - wfx_bom                                                             ) ! Sea-ice mass change through bottom melting 
     177      IF( iom_use('dmtsub') )   CALL iom_put( 'dmtsub', - wfx_sub                                                             ) ! Sea-ice mass change through evaporation and sublimation 
     178      IF( iom_use('dmssub') )   CALL iom_put( 'dmssub', - wfx_snw_sub                                                         ) ! Snow mass change through sublimation 
     179      IF( iom_use('dmisub') )   CALL iom_put( 'dmisub', - wfx_ice_sub                                                         ) ! Sea-ice mass change through sublimation 
     180      IF( iom_use('dmsspr') )   CALL iom_put( 'dmsspr', - wfx_spr                                                             ) ! Snow mass change through snow fall 
     181      IF( iom_use('dmsssi') )   CALL iom_put( 'dmsssi',   wfx_sni*rhos*r1_rhoi                                                ) ! Snow mass change through snow-to-ice conversion 
     182      IF( iom_use('dmsmel') )   CALL iom_put( 'dmsmel', - wfx_snw_sum                                                         ) ! Snow mass change through melt 
     183      IF( iom_use('dmsdyn') )   CALL iom_put( 'dmsdyn', - wfx_snw_dyn + rhos * diag_trp_vs                                    ) ! Snow mass change through dynamics(kg/m2/s) 
     184       
    187185      ! Global ice diagnostics 
    188       IF( iom_use('NH_icearea') .OR. iom_use('NH_icevolu') .OR. iom_use('NH_iceextt') )   THEN   ! NH diagnostics 
    189          ! 
    190          WHERE( ff_t > 0._wp )   ;   zmsk00(:,:) = 1.0e-12 
    191          ELSEWHERE               ;   zmsk00(:,:) = 0. 
    192          END WHERE  
    193          zdiag_area_nh = glob_sum( 'icewri', at_i(:,:) * zmsk00(:,:) * e1e2t(:,:) ) 
    194          zdiag_volu_nh = glob_sum( 'icewri', vt_i(:,:) * zmsk00(:,:) * e1e2t(:,:) ) 
    195          ! 
    196          WHERE( ff_t > 0._wp .AND. at_i > 0.15 )   ; zmsk00(:,:) = 1.0e-12 
    197          ELSEWHERE                                 ; zmsk00(:,:) = 0. 
    198          END WHERE  
    199          zdiag_extt_nh = glob_sum( 'icewri', zmsk00(:,:) * e1e2t(:,:) ) 
    200          ! 
    201          IF( iom_use('NH_icearea') )   CALL iom_put( "NH_icearea" ,  zdiag_area_nh ) 
    202          IF( iom_use('NH_icevolu') )   CALL iom_put( "NH_icevolu" ,  zdiag_volu_nh ) 
    203          IF( iom_use('NH_iceextt') )   CALL iom_put( "NH_iceextt" ,  zdiag_extt_nh ) 
     186      IF(  iom_use('NH_icearea') .OR. iom_use('NH_icevolu') .OR. iom_use('NH_iceextt') .OR. & 
     187         & iom_use('SH_icearea') .OR. iom_use('SH_icevolu') .OR. iom_use('SH_iceextt') ) THEN 
     188         ! 
     189         WHERE( ff_t(:,:) > 0._wp )   ;   z2d(:,:) = 1._wp 
     190         ELSEWHERE                    ;   z2d(:,:) = 0. 
     191         END WHERE 
     192         ! 
     193         IF( iom_use('NH_icearea') )   zdiag_area_nh = glob_sum( 'icewri', at_i *           z2d   * e1e2t * 1.e-12 ) 
     194         IF( iom_use('NH_icevolu') )   zdiag_volu_nh = glob_sum( 'icewri', vt_i *           z2d   * e1e2t * 1.e-12 ) 
     195         IF( iom_use('NH_iceextt') )   zdiag_extt_nh = glob_sum( 'icewri',                  z2d   * e1e2t * 1.e-12 * zmsk15 ) 
     196         ! 
     197         IF( iom_use('SH_icearea') )   zdiag_area_sh = glob_sum( 'icewri', at_i * ( 1._wp - z2d ) * e1e2t * 1.e-12 ) 
     198         IF( iom_use('SH_icevolu') )   zdiag_volu_sh = glob_sum( 'icewri', vt_i * ( 1._wp - z2d ) * e1e2t * 1.e-12 ) 
     199         IF( iom_use('SH_iceextt') )   zdiag_extt_sh = glob_sum( 'icewri',        ( 1._wp - z2d ) * e1e2t * 1.e-12 * zmsk15 ) 
     200         ! 
     201         CALL iom_put( 'NH_icearea' , zdiag_area_nh ) 
     202         CALL iom_put( 'NH_icevolu' , zdiag_volu_nh ) 
     203         CALL iom_put( 'NH_iceextt' , zdiag_extt_nh ) 
     204         CALL iom_put( 'SH_icearea' , zdiag_area_sh ) 
     205         CALL iom_put( 'SH_icevolu' , zdiag_volu_sh ) 
     206         CALL iom_put( 'SH_iceextt' , zdiag_extt_sh ) 
    204207         ! 
    205208      ENDIF 
    206       ! 
    207       IF( iom_use('SH_icearea') .OR. iom_use('SH_icevolu') .OR. iom_use('SH_iceextt') )   THEN   ! SH diagnostics 
    208          ! 
    209          WHERE( ff_t < 0._wp ); zmsk00(:,:) = 1.0e-12;  
    210          ELSEWHERE            ; zmsk00(:,:) = 0. 
    211          END WHERE  
    212          zdiag_area_sh = glob_sum( 'icewri', at_i(:,:) * zmsk00(:,:) * e1e2t(:,:) )  
    213          zdiag_volu_sh = glob_sum( 'icewri', vt_i(:,:) * zmsk00(:,:) * e1e2t(:,:) ) 
    214          ! 
    215          WHERE( ff_t < 0._wp .AND. at_i > 0.15 ); zmsk00(:,:) = 1.0e-12 
    216          ELSEWHERE                              ; zmsk00(:,:) = 0. 
    217          END WHERE  
    218          zdiag_extt_sh = glob_sum( 'icewri', zmsk00(:,:) * e1e2t(:,:) ) 
    219          ! 
    220          IF( iom_use('SH_icearea') ) CALL iom_put( "SH_icearea", zdiag_area_sh ) 
    221          IF( iom_use('SH_icevolu') ) CALL iom_put( "SH_icevolu", zdiag_volu_sh ) 
    222          IF( iom_use('SH_iceextt') ) CALL iom_put( "SH_iceextt", zdiag_extt_sh ) 
    223          ! 
    224       ENDIF  
    225209      ! 
    226210!!CR      !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s 
    227211!!CR      !     IF( kindic < 0 )   CALL ice_wri_state( 'output.abort' ) 
    228212!!CR      !     not yet implemented 
    229 !!gm  idem for the ocean...  Ask Seb how to get read of ioipsl.... 
     213!!gm  idem for the ocean...  Ask Seb how to get rid of ioipsl.... 
    230214      ! 
    231215      IF( ln_timing )  CALL timing_stop('icewri') 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DIA/diadct.F90

    r11380 r11381  
    1111   !!            3.4  ! 09/2011 (C Bricaud) 
    1212   !!---------------------------------------------------------------------- 
    13 #if defined key_diadct 
    14    !!---------------------------------------------------------------------- 
    15    !!   'key_diadct' : 
    16    !!---------------------------------------------------------------------- 
     13   !! 
    1714   !!---------------------------------------------------------------------- 
    1815   !!   dia_dct      :  Compute the transport through a sec. 
     
    4239 
    4340   PUBLIC   dia_dct      ! routine called by step.F90 
    44    PUBLIC   dia_dct_init ! routine called by opa.F90 
    45    PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90  
    46    PRIVATE  readsec 
    47    PRIVATE  removepoints 
    48    PRIVATE  transport 
    49    PRIVATE  dia_dct_wri 
    50  
    51    LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag 
    52  
    53    INTEGER :: nn_dct        ! Frequency of computation 
    54    INTEGER :: nn_dctwri     ! Frequency of output 
    55    INTEGER :: nn_secdebug   ! Number of the section to debug 
     41   PUBLIC   dia_dct_init ! routine called by nemogcm.F90 
     42 
     43   !                         !!** namelist variables ** 
     44   LOGICAL, PUBLIC ::   ln_diadct     ! Calculate transport thru a section or not 
     45   INTEGER         ::   nn_dct        ! Frequency of computation 
     46   INTEGER         ::   nn_dctwri     ! Frequency of output 
     47   INTEGER         ::   nn_secdebug   ! Number of the section to debug 
    5648    
    5749   INTEGER, PARAMETER :: nb_class_max  = 10 
     
    10496CONTAINS 
    10597  
    106   INTEGER FUNCTION diadct_alloc()  
    107      !!----------------------------------------------------------------------  
    108      !!                   ***  FUNCTION diadct_alloc  ***  
    109      !!----------------------------------------------------------------------  
    110      INTEGER :: ierr(2)  
    111      !!----------------------------------------------------------------------  
    112  
    113      ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) )  
    114      ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(2) )  
    115  
    116      diadct_alloc = MAXVAL( ierr )  
    117      IF( diadct_alloc /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' )  
    118   
    119   END FUNCTION diadct_alloc  
    120  
     98   INTEGER FUNCTION diadct_alloc()  
     99      !!----------------------------------------------------------------------  
     100      !!                   ***  FUNCTION diadct_alloc  ***  
     101      !!----------------------------------------------------------------------  
     102 
     103      ALLOCATE( transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), & 
     104         &      transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=diadct_alloc )  
     105 
     106      CALL mpp_sum( 'diadct', diadct_alloc )  
     107      IF( diadct_alloc /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' )  
     108 
     109   END FUNCTION diadct_alloc 
    121110 
    122111   SUBROUTINE dia_dct_init 
     
    130119      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    131120      !! 
    132       NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug 
     121      NAMELIST/nam_diadct/ln_diadct, nn_dct, nn_dctwri, nn_secdebug 
    133122      !!--------------------------------------------------------------------- 
    134123 
    135      REWIND( numnam_ref )              ! Namelist namdct in reference namelist : Diagnostic: transport through sections 
    136      READ  ( numnam_ref, namdct, IOSTAT = ios, ERR = 901) 
    137 901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist' ) 
    138  
    139      REWIND( numnam_cfg )              ! Namelist namdct in configuration namelist : Diagnostic: transport through sections 
    140      READ  ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 ) 
    141 902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist' ) 
    142      IF(lwm) WRITE ( numond, namdct ) 
     124     REWIND( numnam_ref )              ! Namelist nam_diadct in reference namelist : Diagnostic: transport through sections 
     125     READ  ( numnam_ref, nam_diadct, IOSTAT = ios, ERR = 901) 
     126901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadct in reference namelist' ) 
     127 
     128     REWIND( numnam_cfg )              ! Namelist nam_diadct in configuration namelist : Diagnostic: transport through sections 
     129     READ  ( numnam_cfg, nam_diadct, IOSTAT = ios, ERR = 902 ) 
     130902  IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diadct in configuration namelist' ) 
     131     IF(lwm) WRITE ( numond, nam_diadct ) 
    143132 
    144133     IF( lwp ) THEN 
     
    146135        WRITE(numout,*) "diadct_init: compute transports through sections " 
    147136        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" 
    148         WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct 
    149         WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri 
     137        WRITE(numout,*) "       Calculate transport thru sections: ln_diadct = ", ln_diadct 
     138        WRITE(numout,*) "       Frequency of computation:          nn_dct    = ", nn_dct 
     139        WRITE(numout,*) "       Frequency of write:                nn_dctwri = ", nn_dctwri 
    150140 
    151141        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN 
     
    155145        ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug 
    156146        ENDIF 
    157  
     147     ENDIF 
     148 
     149     IF( ln_diadct ) THEN 
     150        ! control 
    158151        IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  & 
    159           &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' ) 
    160  
     152           &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' ) 
     153 
     154        ! allocate dia_dct arrays 
     155        IF( diadct_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' ) 
     156 
     157        !Read section_ijglobal.diadct 
     158        CALL readsec 
     159 
     160        !open output file 
     161        IF( lwm ) THEN 
     162           CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     163           CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     164           CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     165        ENDIF 
     166 
     167        ! Initialise arrays to zero  
     168        transports_3d(:,:,:,:)=0.0  
     169        transports_2d(:,:,:)  =0.0  
     170        ! 
    161171     ENDIF 
    162  
    163      !Read section_ijglobal.diadct 
    164      CALL readsec 
    165  
    166      !open output file 
    167      IF( lwm ) THEN 
    168         CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    169         CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    170         CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    171      ENDIF 
    172  
    173      ! Initialise arrays to zero  
    174      transports_3d(:,:,:,:)=0.0  
    175      transports_2d(:,:,:)  =0.0  
    176172     ! 
    177173  END SUBROUTINE dia_dct_init 
     
    12391235   END FUNCTION interp 
    12401236 
    1241 #else 
    1242    !!---------------------------------------------------------------------- 
    1243    !!   Default option :                                       Dummy module 
    1244    !!---------------------------------------------------------------------- 
    1245    LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag 
    1246    PUBLIC  
    1247    !! $Id$ 
    1248 CONTAINS 
    1249  
    1250    SUBROUTINE dia_dct_init          ! Dummy routine 
    1251       IMPLICIT NONE 
    1252       WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?' 
    1253    END SUBROUTINE dia_dct_init 
    1254  
    1255    SUBROUTINE dia_dct( kt )         ! Dummy routine 
    1256       IMPLICIT NONE 
    1257       INTEGER, INTENT( in ) :: kt   ! ocean time-step index 
    1258       WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt 
    1259    END SUBROUTINE dia_dct 
    1260 #endif 
    1261  
    12621237   !!====================================================================== 
    12631238END MODULE diadct 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DIA/diaharm.F90

    r11380 r11381  
    55   !!====================================================================== 
    66   !! History :  3.1  !  2007  (O. Le Galloudec, J. Chanut)  Original code 
    7    !!---------------------------------------------------------------------- 
    8 #if defined key_diaharm 
    9    !!---------------------------------------------------------------------- 
    10    !!   'key_diaharm' 
    117   !!---------------------------------------------------------------------- 
    128   USE oce             ! ocean dynamics and tracers variables 
     
    2622   IMPLICIT NONE 
    2723   PRIVATE 
    28  
    29    LOGICAL, PUBLIC, PARAMETER :: lk_diaharm  = .TRUE. 
    3024    
    3125   INTEGER, PARAMETER :: jpincomax    = 2.*jpmax_harmo 
     
    3327 
    3428   !                         !!** namelist variables ** 
    35    INTEGER ::   nit000_han    ! First time step used for harmonic analysis 
    36    INTEGER ::   nitend_han    ! Last time step used for harmonic analysis 
    37    INTEGER ::   nstep_han     ! Time step frequency for harmonic analysis 
    38    INTEGER ::   nb_ana        ! Number of harmonics to analyse 
     29   LOGICAL, PUBLIC ::   ln_diaharm    ! Choose tidal harmonic output or not 
     30   INTEGER         ::   nit000_han    ! First time step used for harmonic analysis 
     31   INTEGER         ::   nitend_han    ! Last time step used for harmonic analysis 
     32   INTEGER         ::   nstep_han     ! Time step frequency for harmonic analysis 
     33   INTEGER         ::   nb_ana        ! Number of harmonics to analyse 
    3934 
    4035   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name 
     
    5348   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...) 
    5449 
    55    PUBLIC   dia_harm   ! routine called by step.F90 
     50   PUBLIC   dia_harm        ! routine called by step.F90 
     51   PUBLIC   dia_harm_init   ! routine called by nemogcm.F90 
    5652 
    5753   !!---------------------------------------------------------------------- 
     
    7167      !! 
    7268      !!-------------------------------------------------------------------- 
    73       INTEGER :: jh, nhan, jk, ji 
     69      INTEGER ::   jh, nhan, ji 
    7470      INTEGER ::   ios                 ! Local integer output status for namelist read 
    7571 
    76       NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname 
     72      NAMELIST/nam_diaharm/ ln_diaharm, nit000_han, nitend_han, nstep_han, tname 
    7773      !!---------------------------------------------------------------------- 
    7874 
     
    8278         WRITE(numout,*) '~~~~~~~ ' 
    8379      ENDIF 
    84       ! 
    85       IF( .NOT. ln_tide )   CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis') 
    86       ! 
    87       CALL tide_init_Wave 
    8880      ! 
    8981      REWIND( numnam_ref )              ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis 
     
    9688      ! 
    9789      IF(lwp) THEN 
    98          WRITE(numout,*) 'First time step used for analysis:  nit000_han= ', nit000_han 
    99          WRITE(numout,*) 'Last  time step used for analysis:  nitend_han= ', nitend_han 
    100          WRITE(numout,*) 'Time step frequency for harmonic analysis:  nstep_han= ', nstep_han 
     90         WRITE(numout,*) 'Tidal diagnostics = ', ln_diaharm 
     91         WRITE(numout,*) '   First time step used for analysis:         nit000_han= ', nit000_han 
     92         WRITE(numout,*) '   Last  time step used for analysis:         nitend_han= ', nitend_han 
     93         WRITE(numout,*) '   Time step frequency for harmonic analysis: nstep_han = ', nstep_han 
    10194      ENDIF 
    10295 
    103       ! Basic checks on harmonic analysis time window: 
    104       ! ---------------------------------------------- 
    105       IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   & 
    106          &                                       ' restart capability not implemented' ) 
    107       IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   & 
    108          &                                       'restart capability not implemented' ) 
    109  
    110       IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   & 
    111          &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' ) 
    112  
    113       nb_ana = 0 
    114       DO jk=1,jpmax_harmo 
    115          DO ji=1,jpmax_harmo 
    116             IF(TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN 
    117                nb_ana=nb_ana+1 
    118             ENDIF 
    119          END DO 
    120       END DO 
    121       ! 
    122       IF(lwp) THEN 
    123          WRITE(numout,*) '        Namelist nam_diaharm' 
    124          WRITE(numout,*) '        nb_ana    = ', nb_ana 
    125          CALL flush(numout) 
     96      IF( ln_diaharm .AND. .NOT.ln_tide )   CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis') 
     97 
     98      IF( ln_diaharm ) THEN 
     99 
     100         CALL tide_init_Wave 
     101         ! 
     102         ! Basic checks on harmonic analysis time window: 
     103         ! ---------------------------------------------- 
     104         IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   & 
     105            &                                       ' restart capability not implemented' ) 
     106         IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   & 
     107            &                                       'restart capability not implemented' ) 
     108 
     109         IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   & 
     110            &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' ) 
     111         ! 
     112         nb_ana = 0 
     113         DO jh=1,jpmax_harmo 
     114            DO ji=1,jpmax_harmo 
     115               IF(TRIM(tname(jh)) == Wave(ji)%cname_tide) THEN 
     116                  nb_ana=nb_ana+1 
     117               ENDIF 
     118            END DO 
     119         END DO 
     120         ! 
     121         IF(lwp) THEN 
     122            WRITE(numout,*) '        Namelist nam_diaharm' 
     123            WRITE(numout,*) '        nb_ana    = ', nb_ana 
     124            CALL flush(numout) 
     125         ENDIF 
     126         ! 
     127         IF (nb_ana > jpmax_harmo) THEN 
     128            WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo' 
     129            WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo 
     130            CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 ) 
     131         ENDIF 
     132 
     133         ALLOCATE(name    (nb_ana)) 
     134         DO jh=1,nb_ana 
     135            DO ji=1,jpmax_harmo 
     136               IF (TRIM(tname(jh)) ==  Wave(ji)%cname_tide) THEN 
     137                  name(jh) = ji 
     138                  EXIT 
     139               END IF 
     140            END DO 
     141         END DO 
     142 
     143         ! Initialize frequency array: 
     144         ! --------------------------- 
     145         ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) ) 
     146 
     147         CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana ) 
     148 
     149         IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency ' 
     150 
     151         DO jh = 1, nb_ana 
     152            IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh) 
     153         END DO 
     154 
     155         ! Initialize temporary arrays: 
     156         ! ---------------------------- 
     157         ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 
     158         ana_temp(:,:,:,:) = 0._wp 
     159 
    126160      ENDIF 
    127       ! 
    128       IF (nb_ana > jpmax_harmo) THEN 
    129          WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo' 
    130          WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo 
    131          CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 ) 
    132       ENDIF 
    133  
    134       ALLOCATE(name    (nb_ana)) 
    135       DO jk=1,nb_ana 
    136        DO ji=1,jpmax_harmo 
    137           IF (TRIM(tname(jk)) ==  Wave(ji)%cname_tide) THEN 
    138              name(jk) = ji 
    139              EXIT 
    140           END IF 
    141        END DO 
    142       END DO 
    143  
    144       ! Initialize frequency array: 
    145       ! --------------------------- 
    146       ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) ) 
    147  
    148       CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana ) 
    149  
    150       IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency ' 
    151  
    152       DO jh = 1, nb_ana 
    153         IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh) 
    154       END DO 
    155  
    156       ! Initialize temporary arrays: 
    157       ! ---------------------------- 
    158       ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 
    159       ana_temp(:,:,:,:) = 0._wp 
    160161 
    161162   END SUBROUTINE dia_harm_init 
     
    177178      !!-------------------------------------------------------------------- 
    178179      IF( ln_timing )   CALL timing_start('dia_harm') 
    179       ! 
    180       IF( kt == nit000 )   CALL dia_harm_init 
    181180      ! 
    182181      IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN 
     
    422421      INTEGER, INTENT(in) ::   init  
    423422      ! 
    424       INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd 
     423      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd 
    425424      REAL(wp)                        :: zval1, zval2, zx1 
    426425      REAL(wp), DIMENSION(jpincomax) :: ztmpx, zcol1, zcol2 
     
    434433         ztmp3(:,:) = 0._wp 
    435434         ! 
    436          DO jk1_sd = 1, nsparse 
    437             DO jk2_sd = 1, nsparse 
    438                nisparse(jk2_sd) = nisparse(jk2_sd) 
    439                njsparse(jk2_sd) = njsparse(jk2_sd) 
    440                IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 
    441                   ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  & 
    442                      &                                     + valuesparse(jk1_sd)*valuesparse(jk2_sd) 
     435         DO jh1_sd = 1, nsparse 
     436            DO jh2_sd = 1, nsparse 
     437               nisparse(jh2_sd) = nisparse(jh2_sd) 
     438               njsparse(jh2_sd) = njsparse(jh2_sd) 
     439               IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN 
     440                  ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd))  & 
     441                     &                                     + valuesparse(jh1_sd)*valuesparse(jh2_sd) 
    443442               ENDIF 
    444443            END DO 
     
    515514   END SUBROUTINE SUR_DETERMINE 
    516515 
    517 #else 
    518    !!---------------------------------------------------------------------- 
    519    !!   Default case :   Empty module 
    520    !!---------------------------------------------------------------------- 
    521    LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm = .FALSE. 
    522 CONTAINS 
    523    SUBROUTINE dia_harm ( kt )     ! Empty routine 
    524       INTEGER, INTENT( IN ) :: kt   
    525       WRITE(*,*) 'dia_harm: you should not have seen this print' 
    526    END SUBROUTINE dia_harm 
    527 #endif 
    528  
    529516   !!====================================================================== 
    530517END MODULE diaharm 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DYN/dynspg_ts.F90

    r11380 r11381  
    631631            ! 
    632632         ENDIF     
    633          ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)  
    634          IF ( ln_wd_dl_bc ) THEN 
    635             zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column 
    636             zvwdav2(1:jpi  ,1:jpjm1) = zvwdav2(1:jpi  ,1:jpjm1) + za2 * zvwdmask(1:jpi  ,1:jpjm1)   ! not jpj-row 
    637          END IF 
    638633         !                              
    639634         ! Sum over sub-time-steps to compute advective velocities (only correct on interior domain) 
     
    641636         un_adv(1:jpi,1:jpj) = un_adv(1:jpi,1:jpj) + za2 * zhU(1:jpi,1:jpj) * r1_e2u(1:jpi,1:jpj) 
    642637         vn_adv(1:jpi,1:jpj) = vn_adv(1:jpi,1:jpj) + za2 * zhV(1:jpi,1:jpj) * r1_e1v(1:jpi,1:jpj) 
     638         IF ( ln_wd_dl_bc ) THEN 
     639            zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column 
     640            zvwdav2(1:jpi  ,1:jpjm1) = zvwdav2(1:jpi  ,1:jpjm1) + za2 * zvwdmask(1:jpi  ,1:jpjm1)   ! not jpj-row 
     641         END IF 
    643642         ! 
    644643         ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/IOM/iom.F90

    r11380 r11381  
    5858   PUBLIC iom_init, iom_swap, iom_open, iom_close, iom_setkt, iom_varid, iom_get 
    5959   PUBLIC iom_chkatt, iom_getatt, iom_putatt, iom_getszuld, iom_rstput, iom_delay_rst, iom_put 
    60    PUBLIC iom_use, iom_context_finalize 
     60   PUBLIC iom_use, iom_context_finalize, iom_miss_val 
    6161 
    6262   PRIVATE iom_rp0d, iom_rp1d, iom_rp2d, iom_rp3d 
     
    16711671      CHARACTER(LEN=*), INTENT(in) ::   cdname 
    16721672      REAL(wp)        , INTENT(in) ::   pfield0d 
    1673       REAL(wp)        , DIMENSION(jpi,jpj) ::   zz     ! masson 
     1673!!      REAL(wp)        , DIMENSION(jpi,jpj) ::   zz     ! masson 
    16741674#if defined key_iomput 
    1675       zz(:,:)=pfield0d 
    1676       CALL xios_send_field(cdname, zz) 
    1677       !CALL xios_send_field(cdname, (/pfield0d/))  
     1675!!clem      zz(:,:)=pfield0d 
     1676!!clem      CALL xios_send_field(cdname, zz) 
     1677      CALL xios_send_field(cdname, (/pfield0d/))  
    16781678#else 
    16791679      IF( .FALSE. )   WRITE(numout,*) cdname, pfield0d   ! useless test to avoid compilation warnings 
     
    23912391   !!   NOT 'key_iomput'                               a few dummy routines 
    23922392   !!---------------------------------------------------------------------- 
    2393  
    23942393   SUBROUTINE iom_setkt( kt, cdname ) 
    23952394      INTEGER         , INTENT(in)::   kt  
     
    24062405 
    24072406   LOGICAL FUNCTION iom_use( cdname ) 
    2408       !!---------------------------------------------------------------------- 
    2409       !!---------------------------------------------------------------------- 
    24102407      CHARACTER(LEN=*), INTENT(in) ::   cdname 
    2411       !!---------------------------------------------------------------------- 
    24122408#if defined key_iomput 
    24132409      iom_use = xios_field_is_active( cdname ) 
     
    24162412#endif 
    24172413   END FUNCTION iom_use 
    2418     
     2414 
     2415   SUBROUTINE iom_miss_val( cdname, pmiss_val ) 
     2416      CHARACTER(LEN=*), INTENT(in ) ::   cdname 
     2417      REAL(wp)        , INTENT(out) ::   pmiss_val    
     2418#if defined key_iomput 
     2419      ! get missing value 
     2420      CALL xios_get_field_attr( cdname, default_value = pmiss_val ) 
     2421#else 
     2422      IF( .FALSE. )   WRITE(numout,*) cdname, pmiss_val   ! useless test to avoid compilation warnings 
     2423#endif 
     2424   END SUBROUTINE iom_miss_val 
     2425   
    24192426   !!====================================================================== 
    24202427END MODULE iom 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/nemogcm.F90

    r11380 r11381  
    5959   USE diaobs         ! Observation diagnostics       (dia_obs_init routine) 
    6060   USE diacfl         ! CFL diagnostics               (dia_cfl_init routine) 
     61   USE diaharm        ! tidal harmonics diagnostics  (dia_harm_init routine) 
    6162   USE step           ! NEMO time-stepping                 (stp     routine) 
    6263   USE icbini         ! handle bergs, initialisation 
     
    472473      IF( ln_diacfl    )   CALL dia_cfl_init    ! Initialise CFL diagnostics 
    473474                           CALL dia_ptr_init    ! Poleward TRansports initialization 
    474       IF( lk_diadct    )   CALL dia_dct_init    ! Sections tranports 
     475                           CALL dia_dct_init    ! Sections tranports 
    475476                           CALL dia_hsb_init    ! heat content, salt content and volume budgets 
    476477                           CALL     trd_init    ! Mixed-layer/Vorticity/Integral constraints trends 
     
    478479                           CALL dia_tmb_init    ! TMB outputs 
    479480                           CALL dia_25h_init    ! 25h mean  outputs 
    480       IF( ln_diaobs    )   CALL dia_obs( nit000-1 )   ! Observation operator for restart 
     481                           CALL dia_harm_init   ! tidal harmonics outputs 
     482     IF( ln_diaobs    )    CALL dia_obs( nit000-1 )   ! Observation operator for restart 
    481483 
    482484      !                                      ! Assimilation increments 
     
    639641      USE trc_oce   , ONLY : trc_oce_alloc 
    640642      USE bdy_oce   , ONLY : bdy_oce_alloc 
    641 #if defined key_diadct  
    642       USE diadct    , ONLY : diadct_alloc  
    643 #endif  
    644643      ! 
    645644      INTEGER :: ierr 
     
    653652      ierr = ierr + bdy_oce_alloc()    ! bdy masks (incl. initialization) 
    654653      ! 
    655 #if defined key_diadct  
    656       ierr = ierr + diadct_alloc ()    !  
    657 #endif  
    658       ! 
    659654      CALL mpp_sum( 'nemogcm', ierr ) 
    660655      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'nemo_alloc: unable to allocate standard ocean arrays' ) 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/step.F90

    r11380 r11381  
    217217      IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics 
    218218      IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth) 
    219       IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports 
     219      IF( ln_diadct  )   CALL dia_dct ( kstp )        ! Transports 
    220220                         CALL dia_ar5 ( kstp )        ! ar5 diag 
    221       IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
     221      IF( ln_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
    222222                         CALL dia_wri ( kstp )        ! ocean model: outputs 
    223223      ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/timing.F90

    r11380 r11381  
    657657         ! Compute cpu/elapsed ratio 
    658658         zall_ratio(:) = all_ctime(:) / all_etime(:) 
    659          ztot_ratio    = SUM(zall_ratio(:)) 
    660          zavg_ratio    = ztot_ratio/REAL(jpnij,wp) 
     659         ztot_ratio    = SUM(all_ctime(:))/SUM(all_etime(:)) 
     660         zavg_ratio    = SUM(zall_ratio(:))/REAL(jpnij,wp) 
    661661         zmax_ratio    = MAXVAL(zall_ratio(:)) 
    662662         zmin_ratio    = MINVAL(zall_ratio(:))    
     
    667667         cllignes(2)='1x,"--------------------",//,' 
    668668         cllignes(3)='1x,"Process Rank |"," Elapsed Time (s) |"," CPU Time (s) |"," Ratio CPU/Elapsed",/,' 
    669          cllignes(4)='    (1x,i4,9x,"|",f12.3,6x,"|",f12.3,2x,"|",4x,f7.3,/),' 
    670          WRITE(cllignes(4)(1:4),'(I4)') jpnij 
     669         cllignes(4)='      (4x,i6,4x,"|",f12.3,6x,"|",f12.3,2x,"|",4x,f7.3,/),' 
     670         WRITE(cllignes(4)(1:6),'(I6)') jpnij 
    671671         cllignes(5)='1x,"Total        |",f12.3,6x,"|",F12.3,2x,"|",4x,f7.3,/,' 
    672672         cllignes(6)='1x,"Minimum      |",f12.3,6x,"|",F12.3,2x,"|",4x,f7.3,/,' 
Note: See TracChangeset for help on using the changeset viewer.