Changeset 11381 for NEMO/branches/2019

Timestamp:

2019-07-31T16:44:56+02:00 (5 years 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:

• cfgs/SHARED/namelist_ref (modified) (1 diff)
• src/ICE/ice.F90 (modified) (6 diffs)
• src/ICE/icecor.F90 (modified) (5 diffs)
• src/ICE/icedia.F90 (modified) (4 diffs)
• src/ICE/icedyn.F90 (modified) (4 diffs)
• src/ICE/icedyn_adv.F90 (modified) (1 diff)
• src/ICE/icedyn_rhg.F90 (modified) (1 diff)
• src/ICE/icedyn_rhg_evp.F90 (modified) (20 diffs)
• src/ICE/icesbc.F90 (modified) (3 diffs)
• src/ICE/icestp.F90 (modified) (1 diff)
• src/ICE/iceupdate.F90 (modified) (2 diffs)
• src/ICE/icevar.F90 (modified) (6 diffs)
• src/ICE/icewri.F90 (modified) (4 diffs)
• src/OCE/DIA/diadct.F90 (modified) (7 diffs)
• src/OCE/DIA/diaharm.F90 (modified) (11 diffs)
• src/OCE/DYN/dynspg_ts.F90 (modified) (2 diffs)
• src/OCE/IOM/iom.F90 (modified) (5 diffs)
• src/OCE/nemogcm.F90 (modified) (5 diffs)
• src/OCE/step.F90 (modified) (1 diff)
• src/OCE/timing.F90 (modified) (2 diffs)

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/cfgs/SHARED/namelist_ref

@@ -1306 +1306 @@
    jpni        =   0       !  jpni   number of processors following i (set automatically if < 1)
    jpnj        =   0       !  jpnj   number of processors following j (set automatically if < 1)
+   nn_hlts     =   1       !  added halo width for time splitting
 /
 !-----------------------------------------------------------------------

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/ice.F90

@@ -102 +102 @@
    !! vt_i        |      -      |    Total ice vol. per unit area | m     |
    !! vt_s        |      -      |    Total snow vol. per unit ar. | m     |
+   !! st_i        |      -      |    Total Sea ice salt content   | pss.m |
    !! sm_i        |      -      |    Mean sea ice salinity        | pss   |
    !! tm_i        |      -      |    Mean sea ice temperature     | K     |
@@ -135 +136 @@
    REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice
    LOGICAL , PUBLIC ::   ln_landfast_L16  !: landfast ice parameterizationfrom lemieux2016 
-   LOGICAL , PUBLIC ::   ln_landfast_home !: landfast ice parameterizationfrom home made 
    REAL(wp), PUBLIC ::   rn_depfra        !:    fraction of ocean depth that ice must reach to initiate landfast ice
    REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home) 
@@ -252 +252 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_err_sub !: mass flux error after sublimation                        [kg.m-2.s-1]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total)        [s-1]
-
    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]
    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]
@@ -309 +307 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice !: components of the ice velocity                          (m/s)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i , vt_s  !: ice and snow total volume per unit area                 (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   st_i         !: Total ice salinity content                              (pss.m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i         !: ice total fractional area (ice concentration)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ato_i        !: =1-at_i ; total open water fractional area
@@ -409 +408 @@
          &      wfx_bog    (jpi,jpj) , wfx_dyn   (jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,           &
          &      wfx_res    (jpi,jpj) , wfx_sni   (jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,           &
-         &      afx_tot    (jpi,jpj) , rn_amax_2d(jpi,jpj),                                                  &
+         &      rn_amax_2d (jpi,jpj) ,                                                                       &
          &      qsb_ice_bot(jpi,jpj) , qlead     (jpi,jpj) ,                                                 &
          &      sfx_res    (jpi,jpj) , sfx_bri   (jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  &
@@ -429 +428 @@
       ii = ii + 1
       ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                   &
-         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) ,  &
-         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s (jpi,jpj) ,  &
-         &      sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s (jpi,jpj) ,  &
+         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , st_i(jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) ,  &
+         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s(jpi,jpj) ,  &
+         &      sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s(jpi,jpj) ,  &
          &      om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj)            , STAT=ierr(ii) )
 

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icecor.F90

@@ -84 +84 @@
       !                             !-----------------------------------------------------
       IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        !
-      !                             !-----------------------------------------------------
+         !                          !-----------------------------------------------------
          zzc = rhoi * r1_rdtice
          DO jl = 1, jpl
@@ -117 +117 @@
             END DO
          END DO
-         CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1., v_ice, 'V', -1. )            ! lateral boundary conditions
+         CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1., v_ice, 'V', -1. )
       ENDIF
 
-!!gm I guess the trends are only out on demand 
-!!   So please, only do this is it exite an iom_use of on a these variables
-!!   furthermore, only allocate the diag_ arrays in this case 
-!!   and do the iom_put here so that it is only a local allocation
-!!gm 
       !                             !-----------------------------------------------------
       SELECT CASE( kn )             !  Diagnostics                                       !
@@ -143 +138 @@
          END DO
          !                       ! concentration tendency (dynamics)
-         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 
-         afx_tot(:,:) = zafx(:,:)
-         IF( iom_use('afxdyn') )   CALL iom_put( 'afxdyn' , zafx(:,:) )
+         IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 
+            zafx(:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 
+            CALL iom_put( 'afxdyn' , zafx )
+         ENDIF
          !
       CASE( 2 )                        !--- thermo trend diagnostics & ice aging
@@ -164 +160 @@
          END DO
          !                       ! concentration tendency (total + thermo)
-         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice
-         afx_tot(:,:) = afx_tot(:,:) + zafx(:,:)
-         IF( iom_use('afxthd') )   CALL iom_put( 'afxthd' , zafx(:,:) )
-         IF( iom_use('afxtot') )   CALL iom_put( 'afxtot' , afx_tot(:,:) )
+         IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 
+            zafx(:,:) = zafx(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice
+            CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice )
+            CALL iom_put( 'afxtot' , zafx )
+         ENDIF
          !
       END SELECT
@@ -174 +171 @@
       IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       IF( ln_ctl         )   CALL ice_prt3D   ('icecor')                                                             ! prints
-      IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                   ! prints
+      IF( ln_icectl .AND. kn == 2 ) &
+         &                   CALL ice_prt     ( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                       ! prints
       IF( ln_timing      )   CALL timing_stop ('icecor')                                                             ! timing
       !

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedia.F90

@@ -34 +34 @@
    PUBLIC   ice_dia_init   ! called in icestp.F90
 
-   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents
+   REAL(wp), SAVE ::   z1_e1e2  ! inverse of the ocean area
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini                    ! initial volume, salt and heat contents
    REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends
    
@@ -80 +81 @@
       ENDIF
 
-!!gm glob_sum includes a " * tmask_i ", so remove  " * tmask(:,:,1) "
-
+      IF( kt == nit000 ) THEN
+         z1_e1e2 = 1._wp / glob_sum( 'icedia', e1e2t(:,:) )
+      ENDIF
+      
       ! ----------------------- !
-      ! 1 -  Contents !
+      ! 1 -  Contents           !
       ! ----------------------- !
-      zbg_ivol = glob_sum( 'icedia', vt_i(:,:) * e1e2t(:,:) ) * 1.e-9                  ! ice volume (km3)
-      zbg_svol = glob_sum( 'icedia', vt_s(:,:) * e1e2t(:,:) ) * 1.e-9                  ! snow volume (km3)
-      zbg_area = glob_sum( 'icedia', at_i(:,:) * e1e2t(:,:) ) * 1.e-6                  ! area (km2)
-      zbg_isal = glob_sum( 'icedia', SUM( sv_i(:,:,:), dim=3 ) * e1e2t(:,:) ) * 1.e-9  ! salt content (pss*km3)
-      zbg_item = glob_sum( 'icedia', et_i * e1e2t(:,:) ) * 1.e-20                      ! heat content (1.e20 J)
-      zbg_stem = glob_sum( 'icedia', et_s * e1e2t(:,:) ) * 1.e-20                      ! heat content (1.e20 J)
-      
+      IF(  iom_use('ibgvol_tot' ) .OR. iom_use('sbgvol_tot' ) .OR. iom_use('ibgarea_tot') .OR. &
+         & iom_use('ibgsalt_tot') .OR. iom_use('ibgheat_tot') .OR. iom_use('sbgheat_tot') ) THEN
+
+         zbg_ivol = glob_sum( 'icedia', vt_i(:,:) * e1e2t(:,:) ) * 1.e-9  ! ice volume (km3)
+         zbg_svol = glob_sum( 'icedia', vt_s(:,:) * e1e2t(:,:) ) * 1.e-9  ! snow volume (km3)
+         zbg_area = glob_sum( 'icedia', at_i(:,:) * e1e2t(:,:) ) * 1.e-6  ! area (km2)
+         zbg_isal = glob_sum( 'icedia', st_i(:,:) * e1e2t(:,:) ) * 1.e-9  ! salt content (pss*km3)
+         zbg_item = glob_sum( 'icedia', et_i(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J)
+         zbg_stem = glob_sum( 'icedia', et_s(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J)
+
+         CALL iom_put( 'ibgvol_tot'  , zbg_ivol ) 
+         CALL iom_put( 'sbgvol_tot'  , zbg_svol ) 
+         CALL iom_put( 'ibgarea_tot' , zbg_area ) 
+         CALL iom_put( 'ibgsalt_tot' , zbg_isal ) 
+         CALL iom_put( 'ibgheat_tot' , zbg_item ) 
+         CALL iom_put( 'sbgheat_tot' , zbg_stem ) 
+ 
+      ENDIF
+
       ! ---------------------------!
       ! 2 - Trends due to forcing  !
       ! ---------------------------!
+      ! they must be kept outside an IF(iom_use) because of the call to dia_rst below
       z_frc_volbot = r1_rau0 * glob_sum( 'icedia', -( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-ocean 
       z_frc_voltop = r1_rau0 * glob_sum( 'icedia', -( wfx_sub(:,:) + wfx_spr(:,:) )                    * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-atm
@@ -106 +122 @@
       frc_temtop  = frc_temtop  + z_frc_temtop  * rdt_ice ! 1.e20 J
       frc_tembot  = frc_tembot  + z_frc_tembot  * rdt_ice ! 1.e20 J
+
+      CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water) 
+      CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water) 
+      CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)   
+      CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)   
+      CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)   
+
+      IF(  iom_use('ibgfrchfxtop') .OR. iom_use('ibgfrchfxbot') ) THEN
+         CALL iom_put( 'ibgfrchfxtop' , frc_temtop * z1_e1e2 * 1.e-20 * kt*rdt ) ! heat on top of ice/snw/ocean      (W/m2)
+         CALL iom_put( 'ibgfrchfxbot' , frc_tembot * z1_e1e2 * 1.e-20 * kt*rdt ) ! heat on top of ocean(below ice)   (W/m2) 
+      ENDIF
+      
+      ! ---------------------------------- !
+      ! 3 -  Content variations and drifts !
+      ! ---------------------------------- !
+      IF(  iom_use('ibgvolume') .OR. iom_use('ibgsaltco') .OR. iom_use('ibgheatco') .OR. iom_use('ibgheatfx') ) THEN
             
-      ! ----------------------- !
-      ! 3 -  Content variations !
-      ! ----------------------- !
-      zdiff_vol = r1_rau0 * glob_sum( 'icedia', ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater trend (km3) 
-      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)
-      zdiff_tem =           glob_sum( 'icedia', ( et_i(:,:) + et_s(:,:)           - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20  ! heat content trend (1.e20 J)
-      !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check)
-
-      ! ----------------------- !
-      ! 4 -  Drifts             !
-      ! ----------------------- !
-      zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot )
-      zdiff_sal = zdiff_sal - frc_sal
-      zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop )
-
-      ! ----------------------- !
-      ! 5 - Diagnostics writing !
-      ! ----------------------- !
-!!gm I don't understand the division by the ocean surface (i.e. glob_sum( 'icedia', e1e2t(:,:) ) * 1.e-20 * kt*rdt )
-!!   and its multiplication bu kt ! is it really what we want ? what is this quantity ?
-!!   IF it is really what we want, compute it at kt=nit000, not 3 time by time-step !
-!!   kt*rdt  : you mean rdtice ?
-!!gm
-      !
-      IF( iom_use('ibgvolume')    )   CALL iom_put( 'ibgvolume' , zdiff_vol     )   ! ice/snow volume  drift            (km3 equivalent ocean water)         
-      IF( iom_use('ibgsaltco')    )   CALL iom_put( 'ibgsaltco' , zdiff_sal     )   ! ice salt content drift            (psu*km3 equivalent ocean water)
-      IF( iom_use('ibgheatco')    )   CALL iom_put( 'ibgheatco' , zdiff_tem     )   ! ice/snow heat content drift       (1.e20 J)
-      IF( iom_use('ibgheatfx')    )   CALL iom_put( 'ibgheatfx' ,               &   ! ice/snow heat flux drift          (W/m2)
-         &                                                     zdiff_tem /glob_sum( 'icedia', e1e2t(:,:) * 1.e-20 * kt*rdt ) )
-
-      IF( iom_use('ibgfrcvoltop') )   CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water) 
-      IF( iom_use('ibgfrcvolbot') )   CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water) 
-      IF( iom_use('ibgfrcsal')    )   CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)   
-      IF( iom_use('ibgfrctemtop') )   CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)   
-      IF( iom_use('ibgfrctembot') )   CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)   
-      IF( iom_use('ibgfrchfxtop') )   CALL iom_put( 'ibgfrchfxtop' ,            &   ! heat on top of ice/snw/ocean      (W/m2) 
-         &                                                          frc_temtop / glob_sum( 'icedia', e1e2t(:,:) ) * 1.e-20 * kt*rdt  )
-      IF( iom_use('ibgfrchfxbot') )   CALL iom_put( 'ibgfrchfxbot' ,            &   ! heat on top of ocean(below ice)   (W/m2) 
-         &                                                          frc_tembot / glob_sum( 'icedia', e1e2t(:,:) ) * 1.e-20 * kt*rdt  )
-
-      IF( iom_use('ibgvol_tot' )  )   CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                       (km3)
-      IF( iom_use('sbgvol_tot' )  )   CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                      (km3)
-      IF( iom_use('ibgarea_tot')  )   CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                         (km2)
-      IF( iom_use('ibgsalt_tot')  )   CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content             (pss*km3)
-      IF( iom_use('ibgheat_tot')  )   CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                 (1.e20 J)
-      IF( iom_use('sbgheat_tot')  )   CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                (1.e20 J)
-      !
+         zdiff_vol = r1_rau0 * glob_sum( 'icedia', ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater trend (km3) 
+         zdiff_sal = r1_rau0 * glob_sum( 'icedia', ( rhoi*st_i(:,:)                  - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! salt content trend (km3*pss)
+         zdiff_tem =           glob_sum( 'icedia', ( et_i(:,:) + et_s(:,:)           - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20  ! heat content trend (1.e20 J)
+         !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check)
+         
+         zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot )
+         zdiff_sal = zdiff_sal - frc_sal
+         zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop )
+         
+         CALL iom_put( 'ibgvolume' , zdiff_vol )   ! ice/snow volume  drift            (km3 equivalent ocean water)         
+         CALL iom_put( 'ibgsaltco' , zdiff_sal )   ! ice salt content drift            (psu*km3 equivalent ocean water)
+         CALL iom_put( 'ibgheatco' , zdiff_tem )   ! ice/snow heat content drift       (1.e20 J)
+         !
+      ENDIF
+      
       IF( lrst_ice )   CALL ice_dia_rst( 'WRITE', kt_ice )
       !
@@ -248 +248 @@
             vol_loc_ini(:,:) = rhoi * vt_i(:,:) + rhos * vt_s(:,:)  ! ice/snow volume (kg/m2)
             tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                ! ice/snow heat content (J)
-            sal_loc_ini(:,:) = rhoi * SUM( sv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2)
+            sal_loc_ini(:,:) = rhoi * st_i(:,:)                     ! ice salt content (pss*kg/m2)
          ENDIF
          !

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedyn.F90

@@ -163 +163 @@
             END DO
             CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
-            CALL iom_put( "icediv" , zdivu_i(:,:) )
+            ! output
+            CALL iom_put( 'icediv' , zdivu_i )
+
             DEALLOCATE( zdivu_i )
 
@@ -219 +221 @@
       NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  &
          &             rn_ishlat ,                                                           &
-         &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
+         &             ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
       !!-------------------------------------------------------------------
       !
@@ -242 +244 @@
          WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat
          WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16
-         WRITE(numout,*) '      Landfast: param from home made                         ln_landfast_home= ', ln_landfast_home
          WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra
          WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr
@@ -269 +270 @@
       ENDIF
       !                                      !--- Landfast ice
-      IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp
-      !
-      IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN
-         CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' )
-      ENDIF
+      IF( .NOT.ln_landfast_L16 )   tau_icebfr(:,:) = 0._wp
       !
       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

@@ -100 +100 @@
       diag_trp_vi(:,:) = SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_rdtice
       diag_trp_vs(:,:) = SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_rdtice
-      IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoi          )   ! ice mass transport
-      IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhos          )   ! snw mass transport
-      IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoi * 1.e-03 )   ! salt mass transport (kg/m2/s)
-      IF( iom_use('dihctrp') )   CALL iom_put( "dihctrp" , -diag_trp_ei                )   ! advected ice heat content (W/m2)
-      IF( iom_use('dshctrp') )   CALL iom_put( "dshctrp" , -diag_trp_es                )   ! advected snw heat content (W/m2)
+      IF( iom_use('icemtrp') )   CALL iom_put( 'icemtrp' ,  diag_trp_vi * rhoi          )   ! ice mass transport
+      IF( iom_use('snwmtrp') )   CALL iom_put( 'snwmtrp' ,  diag_trp_vs * rhos          )   ! snw mass transport
+      IF( iom_use('salmtrp') )   CALL iom_put( 'salmtrp' ,  diag_trp_sv * rhoi * 1.e-03 )   ! salt mass transport (kg/m2/s)
+      IF( iom_use('dihctrp') )   CALL iom_put( 'dihctrp' , -diag_trp_ei                 )   ! advected ice heat content (W/m2)
+      IF( iom_use('dshctrp') )   CALL iom_put( 'dshctrp' , -diag_trp_es                 )   ! advected snw heat content (W/m2)
 
       ! controls

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedyn_rhg.F90

@@ -69 +69 @@
          WRITE(numout,*)'ice_dyn_rhg: sea-ice rheology'
          WRITE(numout,*)'~~~~~~~~~~~'
-      ENDIF
-      !
-      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization
-         tau_icebfr(:,:) = 0._wp
-         DO jl = 1, jpl
-            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
-         END DO
       ENDIF
       !

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedyn_rhg_evp.F90

@@ -112 +112 @@
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
       !!
-      LOGICAL, PARAMETER ::   ll_bdy_substep = .TRUE. ! temporary option to call bdy at each sub-time step (T)
-      !                                                                              or only at the main time step (F)
       INTEGER ::   ji, jj       ! dummy loop indices
       INTEGER ::   jter         ! local integers
@@ -137 +135 @@
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ib , ztauV_ib             ! ice-bottom stress at U-V points (landfast param)
-      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
@@ -152 +146 @@
       !                                                                 !    ocean surface (ssh_m) if ice is not embedded
       !                                                                 !    ice bottom surface if ice is embedded   
-      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
-      !
-      REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses
+      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
       REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
 
@@ -163 +162 @@
       REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
       !! --- diags
-      REAL(wp), DIMENSION(jpi,jpj) ::   zswi
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
       !! --- SIMIP diags
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice   
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress  
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
@@ -264 +253 @@
 
       ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
-      IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN   ;   zkt = rn_tensile
-      ELSE                                               ;   zkt = 0._wp
+      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile
+      ELSE                         ;   zkt = 0._wp
       ENDIF
       !
@@ -308 +297 @@
             
             ! Drag ice-atm.
-            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
-            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
+            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
+            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
 
             ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
@@ -316 +305 @@
 
             ! masks
-            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
-            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
+            zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
+            zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
 
             ! switches
-            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zswitchU(ji,jj) = 0._wp
-            ELSE                                                   ;   zswitchU(ji,jj) = 1._wp   ;   ENDIF
-            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zswitchV(ji,jj) = 0._wp
-            ELSE                                                   ;   zswitchV(ji,jj) = 1._wp   ;   ENDIF
+            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
+            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
+            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
+            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
 
          END DO
@@ -339 +328 @@
                ! ice-bottom stress at U points
                zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
-               zTauU_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
+               ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
                ! ice-bottom stress at V points
                zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
-               zTauV_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
+               ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
                ! ice_bottom stress at T points
                zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
-               tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
+               tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
             END DO
          END DO
          CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. )
          !
-      ELSEIF( ln_landfast_home ) THEN          !-- Home made
+      ELSE                               !-- no landfast
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               zTauU_ib(ji,jj) = tau_icebfr(ji,jj)
-               zTauV_ib(ji,jj) = tau_icebfr(ji,jj)
-            END DO
-         END DO
-         !
-      ELSE                                     !-- no landfast
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zTauU_ib(ji,jj) = 0._wp
-               zTauV_ib(ji,jj) = 0._wp
+               ztaux_base(ji,jj) = 0._wp
+               ztauy_base(ji,jj) = 0._wp
             END DO
          END DO
       ENDIF
-      IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) )
 
       !------------------------------------------------------------------------------!
@@ -504 +484 @@
                   !                 !--- tau_bottom/v_ice
                   zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
-                  zTauB = - zTauV_ib(ji,jj) / zvel
+                  zTauB = ztauy_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
                   !
                   !                 !--- Coriolis at V-points (energy conserving formulation)
-                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
                   !
                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+                  zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
-                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
-                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
-                     &             ) * 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
-                     &           ) * zmaskV(ji,jj)
+                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * 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
+                        &           )   * zmsk00y(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
-                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
-                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &              ) * 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
-                     &            ) * zmaskV(ji,jj)
+                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * 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
+                        &            )   * zmsk00y(ji,jj)
                   ENDIF
                END DO
@@ -540 +522 @@
             CALL agrif_interp_ice( 'V' )
 #endif
-            IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' )
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
             !
             DO jj = 2, jpjm1
@@ -552 +534 @@
                   !                 !--- tau_bottom/u_ice
                   zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
-                  zTauB = - zTauU_ib(ji,jj) / zvel
+                  zTauB = ztaux_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
                   !
                   !                 !--- Coriolis at U-points (energy conserving formulation)
-                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
+                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
                   !
                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+                  zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
-                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
-                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
-                     &              ) * 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 
-                     &            ) * zmaskU(ji,jj)
+                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * 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 
+                        &            )   * zmsk00x(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
-                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
-                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &              ) * 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 
-                     &            ) * zmaskU(ji,jj)
+                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * 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 
+                        &            )   * zmsk00x(ji,jj)
                   ENDIF
                END DO
@@ -588 +572 @@
             CALL agrif_interp_ice( 'U' )
 #endif
-            IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' )
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
             !
          ELSE ! odd iterations
@@ -602 +586 @@
                   !                 !--- tau_bottom/u_ice
                   zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
-                  zTauB = - zTauU_ib(ji,jj) / zvel
+                  zTauB = ztaux_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
                   !
                   !                 !--- Coriolis at U-points (energy conserving formulation)
-                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
+                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
                   !
                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+                  zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
-                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
-                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
-                     &              ) * 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 
-                     &            ) * zmaskU(ji,jj)
+                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * 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 
+                        &            )   * zmsk00x(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
-                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
-                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &              ) * 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
-                     &            ) * zmaskU(ji,jj)
+                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * 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
+                        &            )   * zmsk00x(ji,jj)
                   ENDIF
                END DO
@@ -638 +624 @@
             CALL agrif_interp_ice( 'U' )
 #endif
-            IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' )
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
             !
             DO jj = 2, jpjm1
@@ -650 +636 @@
                   !                 !--- tau_bottom/v_ice
                   zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
-                  zTauB = - zTauV_ib(ji,jj) / zvel
+                  zTauB = ztauy_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
                   !
                   !                 !--- Coriolis at v-points (energy conserving formulation)
-                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
                   !
                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+                  zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
-                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
-                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
-                     &             ) * 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
-                     &           ) * zmaskV(ji,jj)
+                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * 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
+                        &           )   * zmsk00y(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
-                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
-                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &              ) * 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
-                     &            ) * zmaskV(ji,jj)
+                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * 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
+                        &            )   * zmsk00y(ji,jj)
                   ENDIF
                END DO
@@ -686 +674 @@
             CALL agrif_interp_ice( 'V' )
 #endif
-            IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' )
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
             !
          ENDIF
@@ -701 +689 @@
       END DO                                              !  end loop over jter  !
       !                                                   ! ==================== !
-      !
-      IF( ln_bdy .AND. .NOT.ll_bdy_substep ) THEN
-         CALL bdy_ice_dyn( 'U' )
-         CALL bdy_ice_dyn( 'V' )
-      ENDIF
       !
       !------------------------------------------------------------------------------!
@@ -764 +747 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
          END DO
       END DO
 
+      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
+      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
+         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
+            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
+         !
+         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
+         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
+         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
+         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
+         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
+         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
+      ENDIF
+       
       ! --- divergence, shear and strength --- !
-      IF( iom_use('icediv') )   CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:) )   ! divergence
-      IF( iom_use('iceshe') )   CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:) )   ! shear
-      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
-
-      ! --- charge ellipse --- !
-      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
+      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
+      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
+      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
+
+      ! --- stress tensor --- !
+      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN
          !
          ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
@@ -780 +778 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
-                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
-                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
+               zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
+                  &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
+                  &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) )
 
                zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress  
 
-               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
+               zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) )
 
 !!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
@@ -799 +797 @@
          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
          !
-         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 )
-         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 )
-         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
-         !
+         CALL iom_put( 'isig1' , zsig1 )
+         CALL iom_put( 'isig2' , zsig2 )
+         CALL iom_put( 'isig3' , zsig3 )
+         !
+         ! Stress tensor invariants (normal and shear stress N/m)
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress
+
          DEALLOCATE( zsig1 , zsig2 , zsig3 )
       ENDIF
       
       ! --- SIMIP --- !
-      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. &
-         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. &
-         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. &
-         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN
-
-         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  &
-            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  &
-            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  &
-            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) )
-         
+      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
+         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., &
+            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )
+
+         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
+         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
+         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
+         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
+         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
+         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
+      ENDIF
+
+      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
+         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
+         !
+         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
+            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
+         !
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
-               
-               ! Stress tensor invariants (normal and shear stress N/m)
-               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
-               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
-               
-               ! Stress terms of the momentum equation (N/m2)
-               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch     ! sea surface slope stress term
-               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
-               
-               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch     ! Coriolis stress term
-               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch
-               
-               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch     ! internal stress term
-               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
-               
-               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
-               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
-               
                ! 2D ice mass, snow mass, area transport arrays (X, Y)
-               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
-               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
-               
+               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
+               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
+
                zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
                zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
-               
+
                zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
                zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
-               
+
                zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
                zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
-               
-            END DO
-         END DO
-         
-         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
-            &                zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
-            &                zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
-            &                zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
-                  
-         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_utau_oi  , 'U', -1., zdiag_vtau_oi  , 'V', -1.,   &
-            &                zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   &
-            &                zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   &
-            &                zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    )
-         
-         IF( iom_use('normstr' ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress
-         IF( iom_use('sheastr' ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress
-         IF( iom_use('dssh_dx' ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x)
-         IF( iom_use('dssh_dy' ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y)
-         IF( iom_use('corstrx' ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x)
-         IF( iom_use('corstry' ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y)
-         IF( iom_use('intstrx' ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x)
-         IF( iom_use('intstry' ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y)
-         IF( iom_use('utau_oi' ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x)
-         IF( iom_use('vtau_oi' ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y)
-         IF( iom_use('xmtrpice') )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s)
-         IF( iom_use('ymtrpice') )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport 
-         IF( iom_use('xmtrpsnw') )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s)
-         IF( iom_use('ymtrpsnw') )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport
-         IF( iom_use('xatrp'   ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport
-         IF( iom_use('yatrp'   ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport
-
-         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  &
-            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  &
-            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  &
-            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     )
+
+            END DO
+         END DO
+
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
+            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
+            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
+
+         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
+         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport 
+         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
+         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
+         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
+         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
+
+         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
+            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
 
       ENDIF

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icesbc.F90

@@ -114 +114 @@
       INTEGER, INTENT(in) ::   ksbc   ! flux formulation (user defined, bulk or Pure Coupled)
       !
-      INTEGER  ::   ji, jj, jl                                ! dummy loop index
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zalb              ! 2D workspace
+      INTEGER  ::   ji, jj, jl      ! dummy loop index
+      REAL(wp) ::   zmiss_val       ! missing value retrieved from xios 
+      REAL(wp), DIMENSION(jpi,jpj,jpl)              ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
+      REAL(wp), DIMENSION(:,:)        , ALLOCATABLE ::   zalb, zmsk00      ! 2D workspace
       !!--------------------------------------------------------------------
       !
@@ -126 +127 @@
          WRITE(numout,*)'~~~~~~~~~~~~~~~'
       ENDIF
+
+      ! get missing value from xml
+      CALL iom_miss_val( "icetemp", zmiss_val )
 
       ! --- cloud-sky and overcast-sky ice albedos --- !
@@ -152 +156 @@
 
       !--- output ice albedo and surface albedo ---!
-      IF( iom_use('icealb') ) THEN
-         WHERE( at_i_b <= epsi06 )   ;   zalb(:,:) = rn_alb_oce
-         ELSEWHERE                   ;   zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
+      IF( iom_use('icealb') .OR. iom_use('albedo') ) THEN
+
+         ALLOCATE( zalb(jpi,jpj), zmsk00(jpi,jpj) )
+
+         WHERE( at_i_b <= epsi06 )
+            zmsk00(:,:) = 0._wp
+            zalb  (:,:) = rn_alb_oce
+         ELSEWHERE
+            zmsk00(:,:) = 1._wp            
+            zalb  (:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
          END WHERE
-         CALL iom_put( "icealb" , zalb(:,:) )
-      ENDIF
-      IF( iom_use('albedo') ) THEN
+         ! ice albedo
+         CALL iom_put( 'icealb' , zalb * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )
+         ! ice+ocean albedo
          zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b )
-         CALL iom_put( "albedo" , zalb(:,:) )
+         CALL iom_put( 'albedo' , zalb )
+
+         DEALLOCATE( zalb, zmsk00 )
+
       ENDIF
       !

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icestp.F90

@@ -425 +425 @@
       wfx_err_sub(:,:) = 0._wp
       !
-      afx_tot(:,:) = 0._wp   ;
-      !
       diag_heat(:,:) = 0._wp ;   diag_sice(:,:) = 0._wp
       diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/iceupdate.F90

@@ -198 +198 @@
       ! --- salt fluxes [kg/m2/s] --- !
       !                           ! sfxice =  sfxbog + sfxbom + sfxsum + sfxsni + sfxopw + sfxres + sfxdyn + sfxbri + sfxsub + sfxlam
-      IF( iom_use('sfxice'  ) )   CALL iom_put( "sfxice", sfx     * 1.e-03 )   ! salt flux from total ice growth/melt
-      IF( iom_use('sfxbog'  ) )   CALL iom_put( "sfxbog", sfx_bog * 1.e-03 )   ! salt flux from bottom growth
-      IF( iom_use('sfxbom'  ) )   CALL iom_put( "sfxbom", sfx_bom * 1.e-03 )   ! salt flux from bottom melting
-      IF( iom_use('sfxsum'  ) )   CALL iom_put( "sfxsum", sfx_sum * 1.e-03 )   ! salt flux from surface melting
-      IF( iom_use('sfxlam'  ) )   CALL iom_put( "sfxlam", sfx_lam * 1.e-03 )   ! salt flux from lateral melting
-      IF( iom_use('sfxsni'  ) )   CALL iom_put( "sfxsni", sfx_sni * 1.e-03 )   ! salt flux from snow ice formation
-      IF( iom_use('sfxopw'  ) )   CALL iom_put( "sfxopw", sfx_opw * 1.e-03 )   ! salt flux from open water formation
-      IF( iom_use('sfxdyn'  ) )   CALL iom_put( "sfxdyn", sfx_dyn * 1.e-03 )   ! salt flux from ridging rafting
-      IF( iom_use('sfxbri'  ) )   CALL iom_put( "sfxbri", sfx_bri * 1.e-03 )   ! salt flux from brines
-      IF( iom_use('sfxres'  ) )   CALL iom_put( "sfxres", sfx_res * 1.e-03 )   ! salt flux from undiagnosed processes
-      IF( iom_use('sfxsub'  ) )   CALL iom_put( "sfxsub", sfx_sub * 1.e-03 )   ! salt flux from sublimation
+      IF( iom_use('sfxice'  ) )   CALL iom_put( 'sfxice', sfx     * 1.e-03 )   ! salt flux from total ice growth/melt
+      IF( iom_use('sfxbog'  ) )   CALL iom_put( 'sfxbog', sfx_bog * 1.e-03 )   ! salt flux from bottom growth
+      IF( iom_use('sfxbom'  ) )   CALL iom_put( 'sfxbom', sfx_bom * 1.e-03 )   ! salt flux from bottom melting
+      IF( iom_use('sfxsum'  ) )   CALL iom_put( 'sfxsum', sfx_sum * 1.e-03 )   ! salt flux from surface melting
+      IF( iom_use('sfxlam'  ) )   CALL iom_put( 'sfxlam', sfx_lam * 1.e-03 )   ! salt flux from lateral melting
+      IF( iom_use('sfxsni'  ) )   CALL iom_put( 'sfxsni', sfx_sni * 1.e-03 )   ! salt flux from snow ice formation
+      IF( iom_use('sfxopw'  ) )   CALL iom_put( 'sfxopw', sfx_opw * 1.e-03 )   ! salt flux from open water formation
+      IF( iom_use('sfxdyn'  ) )   CALL iom_put( 'sfxdyn', sfx_dyn * 1.e-03 )   ! salt flux from ridging rafting
+      IF( iom_use('sfxbri'  ) )   CALL iom_put( 'sfxbri', sfx_bri * 1.e-03 )   ! salt flux from brines
+      IF( iom_use('sfxres'  ) )   CALL iom_put( 'sfxres', sfx_res * 1.e-03 )   ! salt flux from undiagnosed processes
+      IF( iom_use('sfxsub'  ) )   CALL iom_put( 'sfxsub', sfx_sub * 1.e-03 )   ! salt flux from sublimation
 
       ! --- mass fluxes [kg/m2/s] --- !
-      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)
-      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)
+      CALL iom_put( 'emp_oce', emp_oce )   ! emp over ocean (taking into account the snow blown away from the ice)
+      CALL iom_put( 'emp_ice', emp_ice )   ! emp over ice   (taking into account the snow blown away from the ice)
 
       !                           ! vfxice = vfxbog + vfxbom + vfxsum + vfxsni + vfxopw + vfxdyn + vfxres + vfxlam + vfxpnd
-      IF( iom_use('vfxice'  ) )   CALL iom_put( "vfxice" , wfx_ice )   ! mass flux from total ice growth/melt
-      IF( iom_use('vfxbog'  ) )   CALL iom_put( "vfxbog" , wfx_bog )   ! mass flux from bottom growth
-      IF( iom_use('vfxbom'  ) )   CALL iom_put( "vfxbom" , wfx_bom )   ! mass flux from bottom melt 
-      IF( iom_use('vfxsum'  ) )   CALL iom_put( "vfxsum" , wfx_sum )   ! mass flux from surface melt 
-      IF( iom_use('vfxlam'  ) )   CALL iom_put( "vfxlam" , wfx_lam )   ! mass flux from lateral melt 
-      IF( iom_use('vfxsni'  ) )   CALL iom_put( "vfxsni" , wfx_sni )   ! mass flux from snow-ice formation
-      IF( iom_use('vfxopw'  ) )   CALL iom_put( "vfxopw" , wfx_opw )   ! mass flux from growth in open water
-      IF( iom_use('vfxdyn'  ) )   CALL iom_put( "vfxdyn" , wfx_dyn )   ! mass flux from dynamics (ridging)
-      IF( iom_use('vfxres'  ) )   CALL iom_put( "vfxres" , wfx_res )   ! mass flux from undiagnosed processes 
-      IF( iom_use('vfxpnd'  ) )   CALL iom_put( "vfxpnd" , wfx_pnd )   ! mass flux from melt ponds
-      IF( iom_use('vfxsub'  ) )   CALL iom_put( "vfxsub" , wfx_ice_sub )   ! mass flux from ice sublimation (ice-atm.)
-      IF( iom_use('vfxsub_err') ) CALL iom_put( "vfxsub_err", wfx_err_sub )   ! "excess" of sublimation sent to ocean      
-
-      IF ( iom_use( "vfxthin" ) ) THEN   ! mass flux from ice growth in open water + thin ice (<20cm) => comparable to observations  
+      CALL iom_put( 'vfxice'    , wfx_ice     )   ! mass flux from total ice growth/melt
+      CALL iom_put( 'vfxbog'    , wfx_bog     )   ! mass flux from bottom growth
+      CALL iom_put( 'vfxbom'    , wfx_bom     )   ! mass flux from bottom melt 
+      CALL iom_put( 'vfxsum'    , wfx_sum     )   ! mass flux from surface melt 
+      CALL iom_put( 'vfxlam'    , wfx_lam     )   ! mass flux from lateral melt 
+      CALL iom_put( 'vfxsni'    , wfx_sni     )   ! mass flux from snow-ice formation
+      CALL iom_put( 'vfxopw'    , wfx_opw     )   ! mass flux from growth in open water
+      CALL iom_put( 'vfxdyn'    , wfx_dyn     )   ! mass flux from dynamics (ridging)
+      CALL iom_put( 'vfxres'    , wfx_res     )   ! mass flux from undiagnosed processes 
+      CALL iom_put( 'vfxpnd'    , wfx_pnd     )   ! mass flux from melt ponds
+      CALL iom_put( 'vfxsub'    , wfx_ice_sub )   ! mass flux from ice sublimation (ice-atm.)
+      CALL iom_put( 'vfxsub_err', wfx_err_sub )   ! "excess" of sublimation sent to ocean      
+
+      IF ( iom_use( 'vfxthin' ) ) THEN   ! mass flux from ice growth in open water + thin ice (<20cm) => comparable to observations  
          WHERE( hm_i(:,:) < 0.2 .AND. hm_i(:,:) > 0. ) ; z2d = wfx_bog
          ELSEWHERE                                     ; z2d = 0._wp
          END WHERE
-         CALL iom_put( "vfxthin", wfx_opw + z2d )
-      ENDIF
-
-      !                              ! vfxsnw = vfxsnw_sni + vfxsnw_dyn + vfxsnw_sum
-      IF( iom_use('vfxsnw'     ) )   CALL iom_put( "vfxsnw"     , wfx_snw     )   ! mass flux from total snow growth/melt
-      IF( iom_use('vfxsnw_sum' ) )   CALL iom_put( "vfxsnw_sum" , wfx_snw_sum )   ! mass flux from snow melt at the surface
-      IF( iom_use('vfxsnw_sni' ) )   CALL iom_put( "vfxsnw_sni" , wfx_snw_sni )   ! mass flux from snow melt during snow-ice formation 
-      IF( iom_use('vfxsnw_dyn' ) )   CALL iom_put( "vfxsnw_dyn" , wfx_snw_dyn )   ! mass flux from dynamics (ridging) 
-      IF( iom_use('vfxsnw_sub' ) )   CALL iom_put( "vfxsnw_sub" , wfx_snw_sub )   ! mass flux from snow sublimation (ice-atm.) 
-      IF( iom_use('vfxsnw_pre' ) )   CALL iom_put( "vfxsnw_pre" , wfx_spr     )   ! snow precip
+         CALL iom_put( 'vfxthin', wfx_opw + z2d )
+      ENDIF
+
+      !                            ! vfxsnw = vfxsnw_sni + vfxsnw_dyn + vfxsnw_sum
+      CALL iom_put( 'vfxsnw'     , wfx_snw     )   ! mass flux from total snow growth/melt
+      CALL iom_put( 'vfxsnw_sum' , wfx_snw_sum )   ! mass flux from snow melt at the surface
+      CALL iom_put( 'vfxsnw_sni' , wfx_snw_sni )   ! mass flux from snow melt during snow-ice formation 
+      CALL iom_put( 'vfxsnw_dyn' , wfx_snw_dyn )   ! mass flux from dynamics (ridging) 
+      CALL iom_put( 'vfxsnw_sub' , wfx_snw_sub )   ! mass flux from snow sublimation (ice-atm.) 
+      CALL iom_put( 'vfxsnw_pre' , wfx_spr     )   ! snow precip
 
       ! --- heat fluxes [W/m2] --- !
       !                              ! qt_atm_oi - qt_oce_ai = hfxdhc - ( dihctrp + dshctrp )
-      IF( iom_use('qsr_oce'    ) )   CALL iom_put( "qsr_oce"    , qsr_oce * ( 1._wp - at_i_b )                               )   !     solar flux at ocean surface
-      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
-      IF( iom_use('qsr_ice'    ) )   CALL iom_put( "qsr_ice"    , SUM( qsr_ice * a_i_b, dim=3 )                              )   !     solar flux at ice surface
-      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
-      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
-      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
-      IF( iom_use('qt_oce'     ) )   CALL iom_put( "qt_oce"     ,      ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce )
-      IF( iom_use('qt_ice'     ) )   CALL iom_put( "qt_ice"     , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 )     + qemp_ice )
-      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) 
-      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) 
-      IF( iom_use('qemp_oce'   ) )   CALL iom_put( "qemp_oce"   , qemp_oce                                                   )   ! Downward Heat Flux from E-P over ocean
-      IF( iom_use('qemp_ice'   ) )   CALL iom_put( "qemp_ice"   , qemp_ice                                                   )   ! Downward Heat Flux from E-P over ice
+      IF( iom_use('qsr_oce'    ) )   CALL iom_put( 'qsr_oce'    , qsr_oce * ( 1._wp - at_i_b )                               )   !     solar flux at ocean surface
+      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
+      IF( iom_use('qsr_ice'    ) )   CALL iom_put( 'qsr_ice'    , SUM( qsr_ice * a_i_b, dim=3 )                              )   !     solar flux at ice surface
+      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
+      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
+      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
+      IF( iom_use('qt_oce'     ) )   CALL iom_put( 'qt_oce'     ,      ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce )
+      IF( iom_use('qt_ice'     ) )   CALL iom_put( 'qt_ice'     , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 )     + qemp_ice )
+      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) 
+      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) 
+      IF( iom_use('qemp_oce'   ) )   CALL iom_put( 'qemp_oce'   , qemp_oce                                                   )   ! Downward Heat Flux from E-P over ocean
+      IF( iom_use('qemp_ice'   ) )   CALL iom_put( 'qemp_ice'   , qemp_ice                                                   )   ! Downward Heat Flux from E-P over ice
 
       ! heat fluxes from ice transformations
-      !                              ! hfxdhc = hfxbog + hfxbom + hfxsum + hfxopw + hfxdif + hfxsnw - ( hfxthd + hfxdyn + hfxres + hfxsub + hfxspr )
-      IF( iom_use('hfxbog'     ) )   CALL iom_put ("hfxbog"     , hfx_bog             )   ! heat flux used for ice bottom growth 
-      IF( iom_use('hfxbom'     ) )   CALL iom_put ("hfxbom"     , hfx_bom             )   ! heat flux used for ice bottom melt
-      IF( iom_use('hfxsum'     ) )   CALL iom_put ("hfxsum"     , hfx_sum             )   ! heat flux used for ice surface melt
-      IF( iom_use('hfxopw'     ) )   CALL iom_put ("hfxopw"     , hfx_opw             )   ! heat flux used for ice formation in open water
-      IF( iom_use('hfxdif'     ) )   CALL iom_put ("hfxdif"     , hfx_dif             )   ! heat flux used for ice temperature change
-      IF( iom_use('hfxsnw'     ) )   CALL iom_put ("hfxsnw"     , hfx_snw             )   ! heat flux used for snow melt 
-      IF( iom_use('hfxerr'     ) )   CALL iom_put ("hfxerr"     , hfx_err_dif         )   ! heat flux error after heat diffusion (included in qt_oce_ai)
+      !                            ! hfxdhc = hfxbog + hfxbom + hfxsum + hfxopw + hfxdif + hfxsnw - ( hfxthd + hfxdyn + hfxres + hfxsub + hfxspr )
+      CALL iom_put ('hfxbog'     , hfx_bog     )   ! heat flux used for ice bottom growth 
+      CALL iom_put ('hfxbom'     , hfx_bom     )   ! heat flux used for ice bottom melt
+      CALL iom_put ('hfxsum'     , hfx_sum     )   ! heat flux used for ice surface melt
+      CALL iom_put ('hfxopw'     , hfx_opw     )   ! heat flux used for ice formation in open water
+      CALL iom_put ('hfxdif'     , hfx_dif     )   ! heat flux used for ice temperature change
+      CALL iom_put ('hfxsnw'     , hfx_snw     )   ! heat flux used for snow melt 
+      CALL iom_put ('hfxerr'     , hfx_err_dif )   ! heat flux error after heat diffusion (included in qt_oce_ai)
 
       ! heat fluxes associated with mass exchange (freeze/melt/precip...)
-      IF( iom_use('hfxthd'     ) )   CALL iom_put ("hfxthd"     , hfx_thd             )   !  
-      IF( iom_use('hfxdyn'     ) )   CALL iom_put ("hfxdyn"     , hfx_dyn             )   !  
-      IF( iom_use('hfxres'     ) )   CALL iom_put ("hfxres"     , hfx_res             )   !  
-      IF( iom_use('hfxsub'     ) )   CALL iom_put ("hfxsub"     , hfx_sub             )   !  
-      IF( iom_use('hfxspr'     ) )   CALL iom_put ("hfxspr"     , hfx_spr             )   ! Heat flux from snow precip heat content 
+      CALL iom_put ('hfxthd'     , hfx_thd     )   !  
+      CALL iom_put ('hfxdyn'     , hfx_dyn     )   !  
+      CALL iom_put ('hfxres'     , hfx_res     )   !  
+      CALL iom_put ('hfxsub'     , hfx_sub     )   !  
+      CALL iom_put ('hfxspr'     , hfx_spr     )   ! Heat flux from snow precip heat content 
 
       ! other heat fluxes
-      IF( iom_use('hfxsensib'  ) )   CALL iom_put( "hfxsensib"  ,     -qsb_ice_bot * at_i_b         )   ! Sensible oceanic heat flux
-      IF( iom_use('hfxcndbot'  ) )   CALL iom_put( "hfxcndbot"  , SUM( qcn_ice_bot * a_i_b, dim=3 ) )   ! Bottom conduction flux
-      IF( iom_use('hfxcndtop'  ) )   CALL iom_put( "hfxcndtop"  , SUM( qcn_ice_top * a_i_b, dim=3 ) )   ! Surface conduction flux
+      IF( iom_use('hfxsensib'  ) )   CALL iom_put( 'hfxsensib'  ,     -qsb_ice_bot * at_i_b         )   ! Sensible oceanic heat flux
+      IF( iom_use('hfxcndbot'  ) )   CALL iom_put( 'hfxcndbot'  , SUM( qcn_ice_bot * a_i_b, dim=3 ) )   ! Bottom conduction flux
+      IF( iom_use('hfxcndtop'  ) )   CALL iom_put( 'hfxcndtop'  , SUM( qcn_ice_top * a_i_b, dim=3 ) )   ! Surface conduction flux
 
       ! diags
-      IF( iom_use('hfxdhc'     ) )   CALL iom_put ("hfxdhc"     , diag_heat           )   ! Heat content variation in snow and ice 
+      CALL iom_put ('hfxdhc'     , diag_heat   )   ! Heat content variation in snow and ice 
       !
       ! controls
@@ -413 +413 @@
       !! ** Method  :   use of IOM library
       !!----------------------------------------------------------------------
-      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! 'READ'/'WRITE' flag
       INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
       !

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icevar.F90

@@ -32 +32 @@
    !!                        - vt_s(jpi,jpj)
    !!                        - at_i(jpi,jpj)
+   !!                        - st_i(jpi,jpj)
    !!                        - et_s(jpi,jpj)  total snow heat content
    !!                        - et_i(jpi,jpj)  total ice thermal content 
@@ -104 +105 @@
       !
       !                                      ! integrated values
-      vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 )
-      vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 )
-      at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 )
-      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
-      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
+      vt_i(:,:) =       SUM( v_i (:,:,:)           , dim=3 )
+      vt_s(:,:) =       SUM( v_s (:,:,:)           , dim=3 )
+      st_i(:,:) =       SUM( sv_i(:,:,:)           , dim=3 )
+      at_i(:,:) =       SUM( a_i (:,:,:)           , dim=3 )
+      et_s(:,:)  = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 )
+      et_i(:,:)  = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 )
       !
       at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
@@ -138 +140 @@
          tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
          om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
-         sm_i (:,:) = SUM( sv_i(:,:,:)              , dim=3 ) * z1_vt_i(:,:)
+         sm_i (:,:) =      st_i(:,:)                          * z1_vt_i(:,:)
          !
          tm_i(:,:) = 0._wp
@@ -263 +265 @@
       !
       ! integrated values 
-      vt_i (:,:) = SUM( v_i, dim=3 )
-      vt_s (:,:) = SUM( v_s, dim=3 )
-      at_i (:,:) = SUM( a_i, dim=3 )
+      vt_i (:,:) = SUM( v_i , dim=3 )
+      vt_s (:,:) = SUM( v_s , dim=3 )
+      at_i (:,:) = SUM( a_i , dim=3 )
       !
    END SUBROUTINE ice_var_glo2eqv
@@ -533 +535 @@
 
       ! to be sure that at_i is the sum of a_i(jl)
-      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
-      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
+      at_i (:,:) = SUM( a_i (:,:,:), dim=3 )
+      vt_i (:,:) = SUM( v_i (:,:,:), dim=3 )
+!!clem add?
+!      vt_s (:,:) = SUM( v_s (:,:,:), dim=3 )
+!      st_i (:,:) = SUM( sv_i(:,:,:), dim=3 )
+!      et_s(:,:)  = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 )
+!      et_i(:,:)  = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 )
+!!clem
 
       ! open water = 1 if at_i=0
@@ -1085 +1093 @@
          !                              ! ---------------------- !
          CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:) )
-!!$         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:), &
-!!$            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) )
+!!         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:), &
+!!            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) )
          !                              ! ---------------------- !
       ELSEIF( jpl == 1 ) THEN           ! output cat = 1         !
          !                              ! ---------------------- !
-         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1) )         
-!!$         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1), &
-!!$            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) )         
+         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1) )
+!!         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1), &
+!!            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) )         
          !                              ! ----------------------- !
       ELSE                              ! input cat /= output cat !

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icewri.F90

@@ -50 +50 @@
       INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
       REAL(wp) ::   z2da, z2db, zrho1, zrho2
-      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d, zfast !  2D workspace
+      REAL(wp) ::   zmiss_val       ! missing value retrieved from xios 
+      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d, zfast                     ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj)     ::   zmsk00, zmsk05, zmsk15, zmsksn ! O%, 5% and 15% concentration mask and snow mask
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zmsk00l, zmsksnl               ! cat masks
@@ -58 +59 @@
       REAL(wp) ::   zdiag_area_sh, zdiag_extt_sh, zdiag_volu_sh 
       !!-------------------------------------------------------------------
-
+      !
       IF( ln_timing )   CALL timing_start('icewri')
+
+      ! get missing value from xml
+      CALL iom_miss_val( 'icetemp', zmiss_val )
 
       ! brine volume
@@ -85 +89 @@
       ! Standard outputs
       !-----------------
-      zrho1 = ( rau0 - rhoi ) * r1_rau0; zrho2 = rhos * r1_rau0
+      zrho1 = ( rau0 - rhoi ) * r1_rau0 ; zrho2 = rhos * r1_rau0
       ! masks
-      IF( iom_use('icemask'  ) )   CALL iom_put( "icemask"  , zmsk00              )   ! ice mask 0%
-      IF( iom_use('icemask05') )   CALL iom_put( "icemask05", zmsk05              )   ! ice mask 5%
-      IF( iom_use('icemask15') )   CALL iom_put( "icemask15", zmsk15              )   ! ice mask 15%
+      CALL iom_put( 'icemask'  , zmsk00 )   ! ice mask 0%
+      CALL iom_put( 'icemask05', zmsk05 )   ! ice mask 5%
+      CALL iom_put( 'icemask15', zmsk15 )   ! ice mask 15%
+      CALL iom_put( 'icepres'  , zmsk00 )   ! Ice presence (1 or 0) 
       !
       ! general fields
-      IF( iom_use('icemass'  ) )   CALL iom_put( "icemass", rhoi * vt_i * zmsk00  )   ! Ice mass per cell area 
-      IF( iom_use('snwmass'  ) )   CALL iom_put( "snwmass", rhos * vt_s * zmsksn  )   ! Snow mass per cell area
-      IF( iom_use('icepres'  ) )   CALL iom_put( "icepres", zmsk00                )   ! Ice presence (1 or 0) 
-      IF( iom_use('iceconc'  ) )   CALL iom_put( "iceconc", at_i  * zmsk00        )   ! ice concentration
-      IF( iom_use('icevolu'  ) )   CALL iom_put( "icevolu", vt_i  * zmsk00        )   ! ice volume = mean ice thickness over the cell
-      IF( iom_use('icethic'  ) )   CALL iom_put( "icethic", hm_i  * zmsk00        )   ! ice thickness
-      IF( iom_use('snwthic'  ) )   CALL iom_put( "snwthic", hm_s  * zmsk00        )   ! snw thickness
-      IF( iom_use('icebrv'   ) )   CALL iom_put( "icebrv" , bvm_i * zmsk00 * 100. )   ! brine volume
-      IF( iom_use('iceage'   ) )   CALL iom_put( "iceage" , om_i  * zmsk15 / rday )   ! ice age
-      IF( iom_use('icehnew'  ) )   CALL iom_put( "icehnew", ht_i_new              )   ! new ice thickness formed in the leads
-      IF( iom_use('snwvolu'  ) )   CALL iom_put( "snwvolu", vt_s  * zmsksn        )   ! snow volume
-      IF( iom_use('icefrb') ) THEN
+      IF( iom_use('icemass' ) )   CALL iom_put( 'icemass', vt_i * rhoi * zmsk00 )                                           ! Ice mass per cell area 
+      IF( iom_use('snwmass' ) )   CALL iom_put( 'snwmass', vt_s * rhos * zmsksn )                                           ! Snow mass per cell area
+      IF( iom_use('iceconc' ) )   CALL iom_put( 'iceconc', at_i        * zmsk00 )                                           ! ice concentration
+      IF( iom_use('icevolu' ) )   CALL iom_put( 'icevolu', vt_i        * zmsk00 )                                           ! ice volume = mean ice thickness over the cell
+      IF( iom_use('icethic' ) )   CALL iom_put( 'icethic', hm_i        * zmsk00 )                                           ! ice thickness
+      IF( iom_use('snwthic' ) )   CALL iom_put( 'snwthic', hm_s        * zmsk00 )                                           ! snw thickness
+      IF( iom_use('icebrv'  ) )   CALL iom_put( 'icebrv' , bvm_i* 100. * zmsk00 )                                           ! brine volume
+      IF( iom_use('iceage'  ) )   CALL iom_put( 'iceage' , om_i / rday * zmsk15 + zmiss_val * ( 1._wp - zmsk15 ) )          ! ice age
+      IF( iom_use('icehnew' ) )   CALL iom_put( 'icehnew', ht_i_new             )                                           ! new ice thickness formed in the leads
+      IF( iom_use('snwvolu' ) )   CALL iom_put( 'snwvolu', vt_s        * zmsksn )                                           ! snow volume
+      IF( iom_use('icefrb'  ) ) THEN                                                                                        ! Ice freeboard
          z2d(:,:) = ( zrho1 * hm_i(:,:) - zrho2 * hm_s(:,:) )                                         
          WHERE( z2d < 0._wp )   z2d = 0._wp
-                                   CALL iom_put( "icefrb" , z2d * zmsk00          )   ! Ice freeboard
+                                  CALL iom_put( 'icefrb' , z2d * zmsk00         )
       ENDIF
-      !
       ! melt ponds
-      IF( iom_use('iceapnd'  ) )   CALL iom_put( "iceapnd", at_ip  * zmsk00       )   ! melt pond total fraction
-      IF( iom_use('icevpnd'  ) )   CALL iom_put( "icevpnd", vt_ip  * zmsk00       )   ! melt pond total volume per unit area
-      !
+      IF( iom_use('iceapnd' ) )   CALL iom_put( 'iceapnd', at_ip  * zmsk00      )                                           ! melt pond total fraction
+      IF( iom_use('icevpnd' ) )   CALL iom_put( 'icevpnd', vt_ip  * zmsk00      )                                           ! melt pond total volume per unit area
       ! salt
-      IF( iom_use('icesalt'  ) )   CALL iom_put( "icesalt", sm_i  * zmsk00        )   ! mean ice salinity
-      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
-
+      IF( iom_use('icesalt' ) )   CALL iom_put( 'icesalt', sm_i                 * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! mean ice salinity
+      IF( iom_use('icesalm' ) )   CALL iom_put( 'icesalm', st_i * rhoi * 1.0e-3 * zmsk00 )                                  ! Mass of salt in sea ice per cell area
       ! heat
-      IF( iom_use('icetemp'  ) )   CALL iom_put( "icetemp", ( tm_i  - rt0 ) * zmsk00 )   ! ice mean temperature
-      IF( iom_use('snwtemp'  ) )   CALL iom_put( "snwtemp", ( tm_s  - rt0 ) * zmsksn )   ! snw mean temperature
-      IF( iom_use('icettop'  ) )   CALL iom_put( "icettop", ( tm_su - rt0 ) * zmsk00 )   ! temperature at the ice surface
-      IF( iom_use('icetbot'  ) )   CALL iom_put( "icetbot", ( t_bo  - rt0 ) * zmsk00 )   ! temperature at the ice bottom
-      IF( iom_use('icetsni'  ) )   CALL iom_put( "icetsni", ( tm_si - rt0 ) * zmsk00 )   ! temperature at the snow-ice interface
-      IF( iom_use('icehc'    ) )   CALL iom_put( "icehc"  ,  -et_i          * zmsk00 )   ! ice heat content
-      IF( iom_use('snwhc'    ) )   CALL iom_put( "snwhc"  ,  -et_s          * zmsksn )   ! snow heat content
-
+      IF( iom_use('icetemp' ) )   CALL iom_put( 'icetemp', ( tm_i  - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! ice mean temperature
+      IF( iom_use('snwtemp' ) )   CALL iom_put( 'snwtemp', ( tm_s  - rt0 ) * zmsksn + zmiss_val * ( 1._wp - zmsksn ) )      ! snw mean temperature
+      IF( iom_use('icettop' ) )   CALL iom_put( 'icettop', ( tm_su - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the ice surface
+      IF( iom_use('icetbot' ) )   CALL iom_put( 'icetbot', ( t_bo  - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the ice bottom
+      IF( iom_use('icetsni' ) )   CALL iom_put( 'icetsni', ( tm_si - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the snow-ice interface
+      IF( iom_use('icehc'   ) )   CALL iom_put( 'icehc'  ,  -et_i          * zmsk00 )                                       ! ice heat content
+      IF( iom_use('snwhc'   ) )   CALL iom_put( 'snwhc'  ,  -et_s          * zmsksn )                                       ! snow heat content
       ! momentum
-      IF( iom_use('uice'     ) )   CALL iom_put( "uice"   , u_ice                 )   ! ice velocity u component
-      IF( iom_use('vice'     ) )   CALL iom_put( "vice"   , v_ice                 )   ! ice velocity v component
-      IF( iom_use('utau_ai'  ) )   CALL iom_put( "utau_ai", utau_ice * zmsk00     )   ! Wind stress term in force balance (x)
-      IF( iom_use('vtau_ai'  ) )   CALL iom_put( "vtau_ai", vtau_ice * zmsk00     )   ! Wind stress term in force balance (y)
-
-      IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN 
-        ! module of ice velocity
+      IF( iom_use('uice'    ) )   CALL iom_put( 'uice'   , u_ice    )                                                       ! ice velocity u
+      IF( iom_use('vice'    ) )   CALL iom_put( 'vice'   , v_ice    )                                                       ! ice velocity v
+      !
+      IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN                                                              ! module of ice velocity
          DO jj = 2 , jpjm1
             DO ji = 2 , jpim1
-               z2da  = ( u_ice(ji,jj) + u_ice(ji-1,jj) )
-               z2db  = ( v_ice(ji,jj) + v_ice(ji,jj-1) )
+               z2da  = u_ice(ji,jj) + u_ice(ji-1,jj)
+               z2db  = v_ice(ji,jj) + v_ice(ji,jj-1)
                z2d(ji,jj) = 0.5_wp * SQRT( z2da * z2da + z2db * z2db )
            END DO
          END DO
          CALL lbc_lnk( 'icewri', z2d, 'T', 1. )
-         IF( iom_use('icevel') )   CALL iom_put( "icevel" , z2d )
-
-        ! record presence of fast ice
-         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp
+         CALL iom_put( 'icevel', z2d )
+
+         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp                                      ! record presence of fast ice
          ELSEWHERE                                                ; zfast(:,:) = 0._wp
          END WHERE
-         IF( iom_use('fasticepres') )   CALL iom_put( "fasticepres" , zfast )
+         CALL iom_put( 'fasticepres', zfast )
       ENDIF
 
       ! --- category-dependent fields --- !
-      IF( iom_use('icemask_cat' ) )   CALL iom_put( "icemask_cat" , zmsk00l                                                    )   ! ice mask 0%
-      IF( iom_use('iceconc_cat' ) )   CALL iom_put( "iceconc_cat" , a_i * zmsk00l                                              )   ! area for categories
-      IF( iom_use('icethic_cat' ) )   CALL iom_put( "icethic_cat" , h_i * zmsk00l                                              )   ! thickness for categories
-      IF( iom_use('snwthic_cat' ) )   CALL iom_put( "snwthic_cat" , h_s * zmsksnl                                              )   ! snow depth for categories
-      IF( iom_use('icesalt_cat' ) )   CALL iom_put( "icesalt_cat" , s_i * zmsk00l                                              )   ! salinity for categories
-      IF( iom_use('iceage_cat'  ) )   CALL iom_put( "iceage_cat"  , o_i * zmsk00l / rday                                       )   ! ice age
-      IF( iom_use('icetemp_cat' ) )   CALL iom_put( "icetemp_cat" , ( SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i - rt0 ) * zmsk00l )   ! ice temperature
-      IF( iom_use('snwtemp_cat' ) )   CALL iom_put( "snwtemp_cat" , ( SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s - rt0 ) * zmsksnl )   ! snow temperature
-      IF( iom_use('icettop_cat' ) )   CALL iom_put( "icettop_cat" , ( t_su - rt0 ) * zmsk00l                                   )   ! surface temperature
-      IF( iom_use('icebrv_cat'  ) )   CALL iom_put( "icebrv_cat"  ,   bv_i * 100.  * zmsk00l                                   )   ! brine volume
-      IF( iom_use('iceapnd_cat' ) )   CALL iom_put( "iceapnd_cat" ,   a_ip         * zmsk00l                                   )   ! melt pond frac for categories
-      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( "icehpnd_cat" ,   h_ip         * zmsk00l                                   )   ! melt pond frac for categories
-      IF( iom_use('iceafpnd_cat') )   CALL iom_put( "iceafpnd_cat",   a_ip_frac    * zmsk00l                                   )   ! melt pond frac for categories
+      IF( iom_use('icemask_cat' ) )   CALL iom_put( 'icemask_cat' ,                  zmsk00l                                   ) ! ice mask 0%
+      IF( iom_use('iceconc_cat' ) )   CALL iom_put( 'iceconc_cat' , a_i            * zmsk00l                                   ) ! area for categories
+      IF( iom_use('icethic_cat' ) )   CALL iom_put( 'icethic_cat' , h_i            * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! thickness for categories
+      IF( iom_use('snwthic_cat' ) )   CALL iom_put( 'snwthic_cat' , h_s            * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow depth for categories
+      IF( iom_use('icesalt_cat' ) )   CALL iom_put( 'icesalt_cat' , s_i            * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! salinity for categories
+      IF( iom_use('iceage_cat'  ) )   CALL iom_put( 'iceage_cat'  , o_i / rday     * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice age
+      IF( iom_use('icetemp_cat' ) )   CALL iom_put( 'icetemp_cat' , ( SUM( t_i, dim=3 ) * r1_nlay_i - rt0 ) &
+         &                                                                         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice temperature
+      IF( iom_use('snwtemp_cat' ) )   CALL iom_put( 'snwtemp_cat' , ( SUM( t_s, dim=3 ) * r1_nlay_s - rt0 ) &
+         &                                                                         * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow temperature
+      IF( iom_use('icettop_cat' ) )   CALL iom_put( 'icettop_cat' , ( t_su - rt0 ) * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! surface temperature
+      IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume
+      IF( iom_use('iceapnd_cat' ) )   CALL iom_put( 'iceapnd_cat' ,   a_ip         * zmsk00l                                   ) ! melt pond frac for categories
+      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond frac for categories
+      IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories
 
       !------------------
@@ -170 +168 @@
       !------------------
       ! trends
-      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
-      IF( iom_use('dmidyn') )   CALL iom_put( "dmidyn", - wfx_dyn + rhoi * diag_trp_vi      )   ! Sea-ice mass change from dynamics(kg/m2/s)
-      IF( iom_use('dmiopw') )   CALL iom_put( "dmiopw", - wfx_opw                           )   ! Sea-ice mass change through growth in open water
-      IF( iom_use('dmibog') )   CALL iom_put( "dmibog", - wfx_bog                           )   ! Sea-ice mass change through basal growth
-      IF( iom_use('dmisni') )   CALL iom_put( "dmisni", - wfx_sni                           )   ! Sea-ice mass change through snow-to-ice conversion
-      IF( iom_use('dmisum') )   CALL iom_put( "dmisum", - wfx_sum                           )   ! Sea-ice mass change through surface melting
-      IF( iom_use('dmibom') )   CALL iom_put( "dmibom", - wfx_bom                           )   ! Sea-ice mass change through bottom melting
-      IF( iom_use('dmtsub') )   CALL iom_put( "dmtsub", - wfx_sub                           )   ! Sea-ice mass change through evaporation and sublimation
-      IF( iom_use('dmssub') )   CALL iom_put( "dmssub", - wfx_snw_sub                       )   ! Snow mass change through sublimation
-      IF( iom_use('dmisub') )   CALL iom_put( "dmisub", - wfx_ice_sub                       )   ! Sea-ice mass change through sublimation
-      IF( iom_use('dmsspr') )   CALL iom_put( "dmsspr", - wfx_spr                           )   ! Snow mass change through snow fall
-      IF( iom_use('dmsssi') )   CALL iom_put( "dmsssi",   wfx_sni*rhos*r1_rhoi              )   ! Snow mass change through snow-to-ice conversion
-      IF( iom_use('dmsmel') )   CALL iom_put( "dmsmel", - wfx_snw_sum                       )   ! Snow mass change through melt
-      IF( iom_use('dmsdyn') )   CALL iom_put( "dmsdyn", - wfx_snw_dyn + rhos * diag_trp_vs  )   ! Snow mass change through dynamics(kg/m2/s)
-
+      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
+      IF( iom_use('dmidyn') )   CALL iom_put( 'dmidyn', - wfx_dyn + rhoi * diag_trp_vi                                        ) ! Sea-ice mass change from dynamics(kg/m2/s)
+      IF( iom_use('dmiopw') )   CALL iom_put( 'dmiopw', - wfx_opw                                                             ) ! Sea-ice mass change through growth in open water
+      IF( iom_use('dmibog') )   CALL iom_put( 'dmibog', - wfx_bog                                                             ) ! Sea-ice mass change through basal growth
+      IF( iom_use('dmisni') )   CALL iom_put( 'dmisni', - wfx_sni                                                             ) ! Sea-ice mass change through snow-to-ice conversion
+      IF( iom_use('dmisum') )   CALL iom_put( 'dmisum', - wfx_sum                                                             ) ! Sea-ice mass change through surface melting
+      IF( iom_use('dmibom') )   CALL iom_put( 'dmibom', - wfx_bom                                                             ) ! Sea-ice mass change through bottom melting
+      IF( iom_use('dmtsub') )   CALL iom_put( 'dmtsub', - wfx_sub                                                             ) ! Sea-ice mass change through evaporation and sublimation
+      IF( iom_use('dmssub') )   CALL iom_put( 'dmssub', - wfx_snw_sub                                                         ) ! Snow mass change through sublimation
+      IF( iom_use('dmisub') )   CALL iom_put( 'dmisub', - wfx_ice_sub                                                         ) ! Sea-ice mass change through sublimation
+      IF( iom_use('dmsspr') )   CALL iom_put( 'dmsspr', - wfx_spr                                                             ) ! Snow mass change through snow fall
+      IF( iom_use('dmsssi') )   CALL iom_put( 'dmsssi',   wfx_sni*rhos*r1_rhoi                                                ) ! Snow mass change through snow-to-ice conversion
+      IF( iom_use('dmsmel') )   CALL iom_put( 'dmsmel', - wfx_snw_sum                                                         ) ! Snow mass change through melt
+      IF( iom_use('dmsdyn') )   CALL iom_put( 'dmsdyn', - wfx_snw_dyn + rhos * diag_trp_vs                                    ) ! Snow mass change through dynamics(kg/m2/s)
+      
       ! Global ice diagnostics
-      IF( iom_use('NH_icearea') .OR. iom_use('NH_icevolu') .OR. iom_use('NH_iceextt') )   THEN   ! NH diagnostics
-         !
-         WHERE( ff_t > 0._wp )   ;   zmsk00(:,:) = 1.0e-12
-         ELSEWHERE               ;   zmsk00(:,:) = 0.
-         END WHERE 
-         zdiag_area_nh = glob_sum( 'icewri', at_i(:,:) * zmsk00(:,:) * e1e2t(:,:) )
-         zdiag_volu_nh = glob_sum( 'icewri', vt_i(:,:) * zmsk00(:,:) * e1e2t(:,:) )
-         !
-         WHERE( ff_t > 0._wp .AND. at_i > 0.15 )   ; zmsk00(:,:) = 1.0e-12
-         ELSEWHERE                                 ; zmsk00(:,:) = 0.
-         END WHERE 
-         zdiag_extt_nh = glob_sum( 'icewri', zmsk00(:,:) * e1e2t(:,:) )
-         !
-         IF( iom_use('NH_icearea') )   CALL iom_put( "NH_icearea" ,  zdiag_area_nh )
-         IF( iom_use('NH_icevolu') )   CALL iom_put( "NH_icevolu" ,  zdiag_volu_nh )
-         IF( iom_use('NH_iceextt') )   CALL iom_put( "NH_iceextt" ,  zdiag_extt_nh )
+      IF(  iom_use('NH_icearea') .OR. iom_use('NH_icevolu') .OR. iom_use('NH_iceextt') .OR. &
+         & iom_use('SH_icearea') .OR. iom_use('SH_icevolu') .OR. iom_use('SH_iceextt') ) THEN
+         !
+         WHERE( ff_t(:,:) > 0._wp )   ;   z2d(:,:) = 1._wp
+         ELSEWHERE                    ;   z2d(:,:) = 0.
+         END WHERE
+         !
+         IF( iom_use('NH_icearea') )   zdiag_area_nh = glob_sum( 'icewri', at_i *           z2d   * e1e2t * 1.e-12 )
+         IF( iom_use('NH_icevolu') )   zdiag_volu_nh = glob_sum( 'icewri', vt_i *           z2d   * e1e2t * 1.e-12 )
+         IF( iom_use('NH_iceextt') )   zdiag_extt_nh = glob_sum( 'icewri',                  z2d   * e1e2t * 1.e-12 * zmsk15 )
+         !
+         IF( iom_use('SH_icearea') )   zdiag_area_sh = glob_sum( 'icewri', at_i * ( 1._wp - z2d ) * e1e2t * 1.e-12 )
+         IF( iom_use('SH_icevolu') )   zdiag_volu_sh = glob_sum( 'icewri', vt_i * ( 1._wp - z2d ) * e1e2t * 1.e-12 )
+         IF( iom_use('SH_iceextt') )   zdiag_extt_sh = glob_sum( 'icewri',        ( 1._wp - z2d ) * e1e2t * 1.e-12 * zmsk15 )
+         !
+         CALL iom_put( 'NH_icearea' , zdiag_area_nh )
+         CALL iom_put( 'NH_icevolu' , zdiag_volu_nh )
+         CALL iom_put( 'NH_iceextt' , zdiag_extt_nh )
+         CALL iom_put( 'SH_icearea' , zdiag_area_sh )
+         CALL iom_put( 'SH_icevolu' , zdiag_volu_sh )
+         CALL iom_put( 'SH_iceextt' , zdiag_extt_sh )
          !
       ENDIF
-      !
-      IF( iom_use('SH_icearea') .OR. iom_use('SH_icevolu') .OR. iom_use('SH_iceextt') )   THEN   ! SH diagnostics
-         !
-         WHERE( ff_t < 0._wp ); zmsk00(:,:) = 1.0e-12; 
-         ELSEWHERE            ; zmsk00(:,:) = 0.
-         END WHERE 
-         zdiag_area_sh = glob_sum( 'icewri', at_i(:,:) * zmsk00(:,:) * e1e2t(:,:) ) 
-         zdiag_volu_sh = glob_sum( 'icewri', vt_i(:,:) * zmsk00(:,:) * e1e2t(:,:) )
-         !
-         WHERE( ff_t < 0._wp .AND. at_i > 0.15 ); zmsk00(:,:) = 1.0e-12
-         ELSEWHERE                              ; zmsk00(:,:) = 0.
-         END WHERE 
-         zdiag_extt_sh = glob_sum( 'icewri', zmsk00(:,:) * e1e2t(:,:) )
-         !
-         IF( iom_use('SH_icearea') ) CALL iom_put( "SH_icearea", zdiag_area_sh )
-         IF( iom_use('SH_icevolu') ) CALL iom_put( "SH_icevolu", zdiag_volu_sh )
-         IF( iom_use('SH_iceextt') ) CALL iom_put( "SH_iceextt", zdiag_extt_sh )
-         !
-      ENDIF 
       !
 !!CR      !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s
 !!CR      !     IF( kindic < 0 )   CALL ice_wri_state( 'output.abort' )
 !!CR      !     not yet implemented
-!!gm  idem for the ocean...  Ask Seb how to get read of ioipsl....
+!!gm  idem for the ocean...  Ask Seb how to get rid of ioipsl....
       !
       IF( ln_timing )  CALL timing_stop('icewri')

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DIA/diadct.F90

@@ -11 +11 @@
    !!            3.4  ! 09/2011 (C Bricaud)
    !!----------------------------------------------------------------------
-#if defined key_diadct
-   !!----------------------------------------------------------------------
-   !!   'key_diadct' :
-   !!----------------------------------------------------------------------
+   !!
    !!----------------------------------------------------------------------
    !!   dia_dct      :  Compute the transport through a sec.
@@ -42 +39 @@
 
    PUBLIC   dia_dct      ! routine called by step.F90
-   PUBLIC   dia_dct_init ! routine called by opa.F90
-   PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90 
-   PRIVATE  readsec
-   PRIVATE  removepoints
-   PRIVATE  transport
-   PRIVATE  dia_dct_wri
-
-   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
-
-   INTEGER :: nn_dct        ! Frequency of computation
-   INTEGER :: nn_dctwri     ! Frequency of output
-   INTEGER :: nn_secdebug   ! Number of the section to debug
+   PUBLIC   dia_dct_init ! routine called by nemogcm.F90
+
+   !                         !!** namelist variables **
+   LOGICAL, PUBLIC ::   ln_diadct     ! Calculate transport thru a section or not
+   INTEGER         ::   nn_dct        ! Frequency of computation
+   INTEGER         ::   nn_dctwri     ! Frequency of output
+   INTEGER         ::   nn_secdebug   ! Number of the section to debug
    
    INTEGER, PARAMETER :: nb_class_max  = 10
@@ -104 +96 @@
 CONTAINS
  
-  INTEGER FUNCTION diadct_alloc() 
-     !!---------------------------------------------------------------------- 
-     !!                   ***  FUNCTION diadct_alloc  *** 
-     !!---------------------------------------------------------------------- 
-     INTEGER :: ierr(2) 
-     !!---------------------------------------------------------------------- 
-
-     ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) ) 
-     ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(2) ) 
-
-     diadct_alloc = MAXVAL( ierr ) 
-     IF( diadct_alloc /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' ) 
- 
-  END FUNCTION diadct_alloc 
-
+   INTEGER FUNCTION diadct_alloc() 
+      !!---------------------------------------------------------------------- 
+      !!                   ***  FUNCTION diadct_alloc  *** 
+      !!---------------------------------------------------------------------- 
+
+      ALLOCATE( transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), &
+         &      transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=diadct_alloc ) 
+
+      CALL mpp_sum( 'diadct', diadct_alloc ) 
+      IF( diadct_alloc /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' ) 
+
+   END FUNCTION diadct_alloc
 
    SUBROUTINE dia_dct_init
@@ -130 +119 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
+      NAMELIST/nam_diadct/ln_diadct, nn_dct, nn_dctwri, nn_secdebug
       !!---------------------------------------------------------------------
 
-     REWIND( numnam_ref )              ! Namelist namdct in reference namelist : Diagnostic: transport through sections
-     READ  ( numnam_ref, namdct, IOSTAT = ios, ERR = 901)
-901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist' )
-
-     REWIND( numnam_cfg )              ! Namelist namdct in configuration namelist : Diagnostic: transport through sections
-     READ  ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 )
-902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist' )
-     IF(lwm) WRITE ( numond, namdct )
+     REWIND( numnam_ref )              ! Namelist nam_diadct in reference namelist : Diagnostic: transport through sections
+     READ  ( numnam_ref, nam_diadct, IOSTAT = ios, ERR = 901)
+901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadct in reference namelist' )
+
+     REWIND( numnam_cfg )              ! Namelist nam_diadct in configuration namelist : Diagnostic: transport through sections
+     READ  ( numnam_cfg, nam_diadct, IOSTAT = ios, ERR = 902 )
+902  IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diadct in configuration namelist' )
+     IF(lwm) WRITE ( numond, nam_diadct )
 
      IF( lwp ) THEN
@@ -146 +135 @@
         WRITE(numout,*) "diadct_init: compute transports through sections "
         WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
-        WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct
-        WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri
+        WRITE(numout,*) "       Calculate transport thru sections: ln_diadct = ", ln_diadct
+        WRITE(numout,*) "       Frequency of computation:          nn_dct    = ", nn_dct
+        WRITE(numout,*) "       Frequency of write:                nn_dctwri = ", nn_dctwri
 
         IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
@@ -155 +145 @@
         ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug
         ENDIF
-
+     ENDIF
+
+     IF( ln_diadct ) THEN
+        ! control
         IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  &
-          &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
-
+           &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
+
+        ! allocate dia_dct arrays
+        IF( diadct_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' )
+
+        !Read section_ijglobal.diadct
+        CALL readsec
+
+        !open output file
+        IF( lwm ) THEN
+           CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
+           CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
+           CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
+        ENDIF
+
+        ! Initialise arrays to zero 
+        transports_3d(:,:,:,:)=0.0 
+        transports_2d(:,:,:)  =0.0 
+        !
      ENDIF
-
-     !Read section_ijglobal.diadct
-     CALL readsec
-
-     !open output file
-     IF( lwm ) THEN
-        CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
-        CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
-        CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
-     ENDIF
-
-     ! Initialise arrays to zero 
-     transports_3d(:,:,:,:)=0.0 
-     transports_2d(:,:,:)  =0.0 
      !
   END SUBROUTINE dia_dct_init
@@ -1239 +1235 @@
    END FUNCTION interp
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                       Dummy module
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
-   PUBLIC 
-   !! $Id$
-CONTAINS
-
-   SUBROUTINE dia_dct_init          ! Dummy routine
-      IMPLICIT NONE
-      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?'
-   END SUBROUTINE dia_dct_init
-
-   SUBROUTINE dia_dct( kt )         ! Dummy routine
-      IMPLICIT NONE
-      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
-      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
-   END SUBROUTINE dia_dct
-#endif
-
    !!======================================================================
 END MODULE diadct

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DIA/diaharm.F90

@@ -5 +5 @@
    !!======================================================================
    !! History :  3.1  !  2007  (O. Le Galloudec, J. Chanut)  Original code
-   !!----------------------------------------------------------------------
-#if defined key_diaharm
-   !!----------------------------------------------------------------------
-   !!   'key_diaharm'
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
@@ -26 +22 @@
    IMPLICIT NONE
    PRIVATE
-
-   LOGICAL, PUBLIC, PARAMETER :: lk_diaharm  = .TRUE.
    
    INTEGER, PARAMETER :: jpincomax    = 2.*jpmax_harmo
@@ -33 +27 @@
 
    !                         !!** namelist variables **
-   INTEGER ::   nit000_han    ! First time step used for harmonic analysis
-   INTEGER ::   nitend_han    ! Last time step used for harmonic analysis
-   INTEGER ::   nstep_han     ! Time step frequency for harmonic analysis
-   INTEGER ::   nb_ana        ! Number of harmonics to analyse
+   LOGICAL, PUBLIC ::   ln_diaharm    ! Choose tidal harmonic output or not
+   INTEGER         ::   nit000_han    ! First time step used for harmonic analysis
+   INTEGER         ::   nitend_han    ! Last time step used for harmonic analysis
+   INTEGER         ::   nstep_han     ! Time step frequency for harmonic analysis
+   INTEGER         ::   nb_ana        ! Number of harmonics to analyse
 
    INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name
@@ -53 +48 @@
    CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...)
 
-   PUBLIC   dia_harm   ! routine called by step.F90
+   PUBLIC   dia_harm        ! routine called by step.F90
+   PUBLIC   dia_harm_init   ! routine called by nemogcm.F90
 
    !!----------------------------------------------------------------------
@@ -71 +67 @@
       !!
       !!--------------------------------------------------------------------
-      INTEGER :: jh, nhan, jk, ji
+      INTEGER ::   jh, nhan, ji
       INTEGER ::   ios                 ! Local integer output status for namelist read
 
-      NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname
+      NAMELIST/nam_diaharm/ ln_diaharm, nit000_han, nitend_han, nstep_han, tname
       !!----------------------------------------------------------------------
 
@@ -82 +78 @@
          WRITE(numout,*) '~~~~~~~ '
       ENDIF
-      !
-      IF( .NOT. ln_tide )   CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis')
-      !
-      CALL tide_init_Wave
       !
       REWIND( numnam_ref )              ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis
@@ -96 +88 @@
       !
       IF(lwp) THEN
-         WRITE(numout,*) 'First time step used for analysis:  nit000_han= ', nit000_han
-         WRITE(numout,*) 'Last  time step used for analysis:  nitend_han= ', nitend_han
-         WRITE(numout,*) 'Time step frequency for harmonic analysis:  nstep_han= ', nstep_han
+         WRITE(numout,*) 'Tidal diagnostics = ', ln_diaharm
+         WRITE(numout,*) '   First time step used for analysis:         nit000_han= ', nit000_han
+         WRITE(numout,*) '   Last  time step used for analysis:         nitend_han= ', nitend_han
+         WRITE(numout,*) '   Time step frequency for harmonic analysis: nstep_han = ', nstep_han
       ENDIF
 
-      ! Basic checks on harmonic analysis time window:
-      ! ----------------------------------------------
-      IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   &
-         &                                       ' restart capability not implemented' )
-      IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   &
-         &                                       'restart capability not implemented' )
-
-      IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   &
-         &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
-
-      nb_ana = 0
-      DO jk=1,jpmax_harmo
-         DO ji=1,jpmax_harmo
-            IF(TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN
-               nb_ana=nb_ana+1
-            ENDIF
-         END DO
-      END DO
-      !
-      IF(lwp) THEN
-         WRITE(numout,*) '        Namelist nam_diaharm'
-         WRITE(numout,*) '        nb_ana    = ', nb_ana
-         CALL flush(numout)
+      IF( ln_diaharm .AND. .NOT.ln_tide )   CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis')
+
+      IF( ln_diaharm ) THEN
+
+         CALL tide_init_Wave
+         !
+         ! Basic checks on harmonic analysis time window:
+         ! ----------------------------------------------
+         IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   &
+            &                                       ' restart capability not implemented' )
+         IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   &
+            &                                       'restart capability not implemented' )
+
+         IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   &
+            &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
+         !
+         nb_ana = 0
+         DO jh=1,jpmax_harmo
+            DO ji=1,jpmax_harmo
+               IF(TRIM(tname(jh)) == Wave(ji)%cname_tide) THEN
+                  nb_ana=nb_ana+1
+               ENDIF
+            END DO
+         END DO
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) '        Namelist nam_diaharm'
+            WRITE(numout,*) '        nb_ana    = ', nb_ana
+            CALL flush(numout)
+         ENDIF
+         !
+         IF (nb_ana > jpmax_harmo) THEN
+            WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo'
+            WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo
+            CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 )
+         ENDIF
+
+         ALLOCATE(name    (nb_ana))
+         DO jh=1,nb_ana
+            DO ji=1,jpmax_harmo
+               IF (TRIM(tname(jh)) ==  Wave(ji)%cname_tide) THEN
+                  name(jh) = ji
+                  EXIT
+               END IF
+            END DO
+         END DO
+
+         ! Initialize frequency array:
+         ! ---------------------------
+         ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
+
+         CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
+
+         IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
+
+         DO jh = 1, nb_ana
+            IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh)
+         END DO
+
+         ! Initialize temporary arrays:
+         ! ----------------------------
+         ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
+         ana_temp(:,:,:,:) = 0._wp
+
       ENDIF
-      !
-      IF (nb_ana > jpmax_harmo) THEN
-         WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo'
-         WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo
-         CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 )
-      ENDIF
-
-      ALLOCATE(name    (nb_ana))
-      DO jk=1,nb_ana
-       DO ji=1,jpmax_harmo
-          IF (TRIM(tname(jk)) ==  Wave(ji)%cname_tide) THEN
-             name(jk) = ji
-             EXIT
-          END IF
-       END DO
-      END DO
-
-      ! Initialize frequency array:
-      ! ---------------------------
-      ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
-
-      CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
-
-      IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
-
-      DO jh = 1, nb_ana
-        IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh)
-      END DO
-
-      ! Initialize temporary arrays:
-      ! ----------------------------
-      ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
-      ana_temp(:,:,:,:) = 0._wp
 
    END SUBROUTINE dia_harm_init
@@ -177 +178 @@
       !!--------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('dia_harm')
-      !
-      IF( kt == nit000 )   CALL dia_harm_init
       !
       IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
@@ -422 +421 @@
       INTEGER, INTENT(in) ::   init 
       !
-      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
+      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd
       REAL(wp)                        :: zval1, zval2, zx1
       REAL(wp), DIMENSION(jpincomax) :: ztmpx, zcol1, zcol2
@@ -434 +433 @@
          ztmp3(:,:) = 0._wp
          !
-         DO jk1_sd = 1, nsparse
-            DO jk2_sd = 1, nsparse
-               nisparse(jk2_sd) = nisparse(jk2_sd)
-               njsparse(jk2_sd) = njsparse(jk2_sd)
-               IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
-                  ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  &
-                     &                                     + valuesparse(jk1_sd)*valuesparse(jk2_sd)
+         DO jh1_sd = 1, nsparse
+            DO jh2_sd = 1, nsparse
+               nisparse(jh2_sd) = nisparse(jh2_sd)
+               njsparse(jh2_sd) = njsparse(jh2_sd)
+               IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN
+                  ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd))  &
+                     &                                     + valuesparse(jh1_sd)*valuesparse(jh2_sd)
                ENDIF
             END DO
@@ -515 +514 @@
    END SUBROUTINE SUR_DETERMINE
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Default case :   Empty module
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm = .FALSE.
-CONTAINS
-   SUBROUTINE dia_harm ( kt )     ! Empty routine
-      INTEGER, INTENT( IN ) :: kt  
-      WRITE(*,*) 'dia_harm: you should not have seen this print'
-   END SUBROUTINE dia_harm
-#endif
-
    !!======================================================================
 END MODULE diaharm

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DYN/dynspg_ts.F90

@@ -631 +631 @@
             !
          ENDIF    
-         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 
-         IF ( ln_wd_dl_bc ) THEN
-            zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column
-            zvwdav2(1:jpi  ,1:jpjm1) = zvwdav2(1:jpi  ,1:jpjm1) + za2 * zvwdmask(1:jpi  ,1:jpjm1)   ! not jpj-row
-         END IF
          !                             
          ! Sum over sub-time-steps to compute advective velocities (only correct on interior domain)
@@ -641 +636 @@
          un_adv(1:jpi,1:jpj) = un_adv(1:jpi,1:jpj) + za2 * zhU(1:jpi,1:jpj) * r1_e2u(1:jpi,1:jpj)
          vn_adv(1:jpi,1:jpj) = vn_adv(1:jpi,1:jpj) + za2 * zhV(1:jpi,1:jpj) * r1_e1v(1:jpi,1:jpj)
+         IF ( ln_wd_dl_bc ) THEN
+            zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column
+            zvwdav2(1:jpi  ,1:jpjm1) = zvwdav2(1:jpi  ,1:jpjm1) + za2 * zvwdmask(1:jpi  ,1:jpjm1)   ! not jpj-row
+         END IF
          !
          !

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/IOM/iom.F90

@@ -58 +58 @@
    PUBLIC iom_init, iom_swap, iom_open, iom_close, iom_setkt, iom_varid, iom_get
    PUBLIC iom_chkatt, iom_getatt, iom_putatt, iom_getszuld, iom_rstput, iom_delay_rst, iom_put
-   PUBLIC iom_use, iom_context_finalize
+   PUBLIC iom_use, iom_context_finalize, iom_miss_val
 
    PRIVATE iom_rp0d, iom_rp1d, iom_rp2d, iom_rp3d
@@ -1671 +1671 @@
       CHARACTER(LEN=*), INTENT(in) ::   cdname
       REAL(wp)        , INTENT(in) ::   pfield0d
-      REAL(wp)        , DIMENSION(jpi,jpj) ::   zz     ! masson
+!!      REAL(wp)        , DIMENSION(jpi,jpj) ::   zz     ! masson
 #if defined key_iomput
-      zz(:,:)=pfield0d
-      CALL xios_send_field(cdname, zz)
-      !CALL xios_send_field(cdname, (/pfield0d/)) 
+!!clem      zz(:,:)=pfield0d
+!!clem      CALL xios_send_field(cdname, zz)
+      CALL xios_send_field(cdname, (/pfield0d/)) 
 #else
       IF( .FALSE. )   WRITE(numout,*) cdname, pfield0d   ! useless test to avoid compilation warnings
@@ -2391 +2391 @@
    !!   NOT 'key_iomput'                               a few dummy routines
    !!----------------------------------------------------------------------
-
    SUBROUTINE iom_setkt( kt, cdname )
       INTEGER         , INTENT(in)::   kt 
@@ -2406 +2405 @@
 
    LOGICAL FUNCTION iom_use( cdname )
-      !!----------------------------------------------------------------------
-      !!----------------------------------------------------------------------
       CHARACTER(LEN=*), INTENT(in) ::   cdname
-      !!----------------------------------------------------------------------
 #if defined key_iomput
       iom_use = xios_field_is_active( cdname )
@@ -2416 +2412 @@
 #endif
    END FUNCTION iom_use
-   
+
+   SUBROUTINE iom_miss_val( cdname, pmiss_val )
+      CHARACTER(LEN=*), INTENT(in ) ::   cdname
+      REAL(wp)        , INTENT(out) ::   pmiss_val   
+#if defined key_iomput
+      ! get missing value
+      CALL xios_get_field_attr( cdname, default_value = pmiss_val )
+#else
+      IF( .FALSE. )   WRITE(numout,*) cdname, pmiss_val   ! useless test to avoid compilation warnings
+#endif
+   END SUBROUTINE iom_miss_val
+  
    !!======================================================================
 END MODULE iom

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/nemogcm.F90

@@ -59 +59 @@
    USE diaobs         ! Observation diagnostics       (dia_obs_init routine)
    USE diacfl         ! CFL diagnostics               (dia_cfl_init routine)
+   USE diaharm        ! tidal harmonics diagnostics  (dia_harm_init routine)
    USE step           ! NEMO time-stepping                 (stp     routine)
    USE icbini         ! handle bergs, initialisation
@@ -472 +473 @@
       IF( ln_diacfl    )   CALL dia_cfl_init    ! Initialise CFL diagnostics
                            CALL dia_ptr_init    ! Poleward TRansports initialization
-      IF( lk_diadct    )   CALL dia_dct_init    ! Sections tranports
+                           CALL dia_dct_init    ! Sections tranports
                            CALL dia_hsb_init    ! heat content, salt content and volume budgets
                            CALL     trd_init    ! Mixed-layer/Vorticity/Integral constraints trends
@@ -478 +479 @@
                            CALL dia_tmb_init    ! TMB outputs
                            CALL dia_25h_init    ! 25h mean  outputs
-      IF( ln_diaobs    )   CALL dia_obs( nit000-1 )   ! Observation operator for restart
+                           CALL dia_harm_init   ! tidal harmonics outputs
+     IF( ln_diaobs    )    CALL dia_obs( nit000-1 )   ! Observation operator for restart
 
       !                                      ! Assimilation increments
@@ -639 +641 @@
       USE trc_oce   , ONLY : trc_oce_alloc
       USE bdy_oce   , ONLY : bdy_oce_alloc
-#if defined key_diadct 
-      USE diadct    , ONLY : diadct_alloc 
-#endif 
       !
       INTEGER :: ierr
@@ -653 +652 @@
       ierr = ierr + bdy_oce_alloc()    ! bdy masks (incl. initialization)
       !
-#if defined key_diadct 
-      ierr = ierr + diadct_alloc ()    ! 
-#endif 
-      !
       CALL mpp_sum( 'nemogcm', ierr )
       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

@@ -217 +217 @@
       IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics
       IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth)
-      IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports
+      IF( ln_diadct  )   CALL dia_dct ( kstp )        ! Transports
                          CALL dia_ar5 ( kstp )        ! ar5 diag
-      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
+      IF( ln_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
                          CALL dia_wri ( kstp )        ! ocean model: outputs
       !

NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/timing.F90

@@ -657 +657 @@
          ! Compute cpu/elapsed ratio
          zall_ratio(:) = all_ctime(:) / all_etime(:)
-         ztot_ratio    = SUM(zall_ratio(:))
-         zavg_ratio    = ztot_ratio/REAL(jpnij,wp)
+         ztot_ratio    = SUM(all_ctime(:))/SUM(all_etime(:))
+         zavg_ratio    = SUM(zall_ratio(:))/REAL(jpnij,wp)
          zmax_ratio    = MAXVAL(zall_ratio(:))
          zmin_ratio    = MINVAL(zall_ratio(:))   
@@ -667 +667 @@
          cllignes(2)='1x,"--------------------",//,'
          cllignes(3)='1x,"Process Rank |"," Elapsed Time (s) |"," CPU Time (s) |"," Ratio CPU/Elapsed",/,'
-         cllignes(4)='    (1x,i4,9x,"|",f12.3,6x,"|",f12.3,2x,"|",4x,f7.3,/),'
-         WRITE(cllignes(4)(1:4),'(I4)') jpnij
+         cllignes(4)='      (4x,i6,4x,"|",f12.3,6x,"|",f12.3,2x,"|",4x,f7.3,/),'
+         WRITE(cllignes(4)(1:6),'(I6)') jpnij
          cllignes(5)='1x,"Total        |",f12.3,6x,"|",F12.3,2x,"|",4x,f7.3,/,'
          cllignes(6)='1x,"Minimum      |",f12.3,6x,"|",F12.3,2x,"|",4x,f7.3,/,'