Ignore:
Timestamp:
2017-09-12T20:46:13+02:00 (3 years ago)
Author:
clem
Message:

changes in style - part6 - one more round

File:
1 edited

Legend:

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

    r8514 r8517  
    4747CONTAINS 
    4848 
    49    SUBROUTINE ice_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 
     49   SUBROUTINE ice_cons_hsm( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft ) 
    5050      !!---------------------------------------------------------------------- 
    5151      !!                       ***  ROUTINE ice_cons_hsm *** 
     
    5656      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
    5757      !!              It prints in ocean.output if there is a violation of conservation at each time-step 
    58       !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to 
     58      !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine violations are set to 
    5959      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
    6060      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
     
    6363      INTEGER         , INTENT(in)    ::   icount        ! called at: =0 the begining of the routine, =1  the end 
    6464      CHARACTER(len=*), INTENT(in)    ::   cd_routine    ! name of the routine 
    65       REAL(wp)        , INTENT(inout) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b   ! ???? 
    66       !! 
    67       REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,  zft 
    68       REAL(wp)                        :: zvmin, zamin, zamax  
    69       REAL(wp)                        :: zvtrp, zetrp 
    70       REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill 
    71       REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
     65      REAL(wp)        , INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 
     66      !! 
     67      REAL(wp) ::   zv, zs, zt, zfs, zfv, zft 
     68      REAL(wp) ::  zvmin, zamin, zamax  
     69      REAL(wp) ::  zvtrp, zetrp 
     70      REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill 
     71      REAL(wp), PARAMETER ::  zconv = 1.e-9 ! convert W to GW and kg to Mt 
    7272      !!---------------------------------------------------------------------- 
    7373      ! 
    7474      IF( icount == 0 ) THEN 
     75         !                          ! water flux 
     76         pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
     77            &                    wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
     78            &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    & 
     79            &                  ) * e1e2t(:,:) ) * zconv 
     80         ! 
    7581         !                          ! salt flux 
    76          zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    77             &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
    78             &                ) *  e1e2t(:,:) ) * zconv  
    79          ! 
    80          !                          ! water flux 
    81          zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
    82             &                  wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
    83             &                  wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    & 
    84             &                ) * e1e2t(:,:) ) * zconv 
     82         pdiag_fs = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
     83            &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
     84            &                  ) *  e1e2t(:,:) ) * zconv  
    8585         ! 
    8686         !                          ! heat flux 
    87          zft_b = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    88             &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    89             &                ) *  e1e2t(:,:) ) * zconv 
    90  
    91          zvi_b = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv ) 
    92  
    93          zsmv_b = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t * zconv ) 
    94  
    95          zei_b = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
    96             &                + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv 
     87         pdiag_ft = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
     88            &                  - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
     89            &                  ) *  e1e2t(:,:) ) * zconv 
     90 
     91         pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv ) 
     92 
     93         pdiag_s = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t * zconv ) 
     94 
     95         pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
     96            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv 
    9797 
    9898      ELSEIF( icount == 1 ) THEN 
     99 
     100         ! water flux 
     101         zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
     102            &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
     103            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        & 
     104            &              ) * e1e2t(:,:) ) * zconv - pdiag_fv 
    99105 
    100106         ! salt flux 
    101107         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    102108            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &  
    103             &              ) * e1e2t(:,:) ) * zconv - zfs_b 
    104  
    105          ! water flux 
    106          zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
    107             &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
    108             &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        & 
    109             &              ) * e1e2t(:,:) ) * zconv - zfw_b 
     109            &              ) * e1e2t(:,:) ) * zconv - pdiag_fs 
    110110 
    111111         ! heat flux 
    112112         zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    113113            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    114             &              ) * e1e2t(:,:) ) * zconv - zft_b 
     114            &              ) * e1e2t(:,:) ) * zconv - pdiag_ft 
    115115  
    116116         ! outputs 
    117          zvi = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  & 
    118             &       - zvi_b ) * r1_rdtice - zfw ) * rday 
    119  
    120          zsmv = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  & 
    121             &       - zsmv_b ) * r1_rdtice + zfs ) * rday 
    122  
    123          zei = ( glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )   & 
     117         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  & 
     118            &     - pdiag_v ) * r1_rdtice - zfv ) * rday 
     119 
     120         zs = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  & 
     121            &     - pdiag_s ) * r1_rdtice + zfs ) * rday 
     122 
     123         zt = ( glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )   & 
    124124            &                + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv   & 
    125             &   - zei_b ) * r1_rdtice + zft 
     125            &   - pdiag_t ) * r1_rdtice + zft 
    126126 
    127127         ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 
     
    137137         zv_sill = zarea * 2.5e-5 
    138138         zs_sill = zarea * 25.e-5 
    139          zh_sill = zarea * 10.e-5 
     139         zt_sill = zarea * 10.e-5 
    140140 
    141141         IF(lwp) THEN 
    142             IF ( ABS( zvi  ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zvi 
    143             IF ( ABS( zsmv ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zsmv 
    144             IF ( ABS( zei  ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zei 
     142            IF ( ABS( zv   ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv 
     143            IF ( ABS( zs   ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 
     144            IF ( ABS( zt   ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt 
    145145            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'iceadv' ) THEN 
    146                                          WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp 
    147                                          WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp 
    148             ENDIF 
    149             IF (     zvmin   < -epsi10 ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
    150             IF (     zamax   > MAX( rn_amax_n, rn_amax_s ) + epsi10 .AND. & 
    151                &                         cd_routine /= 'iceadv' .AND. cd_routine /= 'icerdgrft' ) THEN 
    152                                          WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
    153             IF (     zamax   > 1._wp   ) WRITE(numout,*) 'violation a_i>1               (',cd_routine,') = ',zamax 
    154             ENDIF 
    155             IF (      zamin  < -epsi10 ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
     146                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp 
     147                                           WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp 
     148            ENDIF 
     149            IF ( zvmin < -epsi10 )         WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
     150            IF ( zamax > MAX(rn_amax_n,rn_amax_s) + epsi10 .AND. cd_routine /= 'iceadv' .AND. cd_routine /= 'icerdgrft' )  & 
     151               &                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
     152            IF ( zamin < -epsi10 )         WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    156153         ENDIF 
    157154         ! 
     
    169166      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
    170167      !!              It prints in ocean.output if there is a violation of conservation at each time-step 
    171       !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to 
     168      !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine the violation are set to 
    172169      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
    173170      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
     
    176173      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine 
    177174      REAL(wp)                        :: zhfx, zsfx, zvfx 
    178       REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill 
     175      REAL(wp)                        :: zarea, zv_sill, zs_sill, zt_sill 
    179176      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
    180177      !!---------------------------------------------------------------------- 
     178 
     179      ! water flux 
     180      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday 
     181 
     182      ! salt flux 
     183      zsfx  = glob_sum( ( sfx + diag_smvi ) * e1e2t ) * zconv * rday 
    181184 
    182185      ! heat flux 
     
    184187      !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check) 
    185188         &              ) * e1e2t ) * zconv 
    186       ! salt flux 
    187       zsfx  = glob_sum( ( sfx + diag_smvi ) * e1e2t ) * zconv * rday 
    188       ! water flux 
    189       zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday 
    190189 
    191190      ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     
    193192      zv_sill = zarea * 2.5e-5 
    194193      zs_sill = zarea * 25.e-5 
    195       zh_sill = zarea * 10.e-5 
     194      zt_sill = zarea * 10.e-5 
    196195 
    197196      IF(lwp) THEN 
    198          IF( ABS( zvfx ) > zv_sill )   WRITE(numout,*) 'violation vfx    [Mt/day]       (',cd_routine,')  = ',(zvfx) 
    199          IF( ABS( zsfx ) > zs_sill )   WRITE(numout,*) 'violation sfx    [psu*Mt/day]   (',cd_routine,')  = ',(zsfx) 
    200          IF( ABS( zhfx ) > zh_sill )   WRITE(numout,*) 'violation hfx    [GW]           (',cd_routine,')  = ',(zhfx) 
     197         IF( ABS( zvfx ) > zv_sill )   WRITE(numout,*) 'violation vfx  [Mt/day]       (',cd_routine,') = ',zvfx 
     198         IF( ABS( zsfx ) > zs_sill )   WRITE(numout,*) 'violation sfx  [psu*Mt/day]   (',cd_routine,') = ',zsfx 
     199         IF( ABS( zhfx ) > zt_sill )   WRITE(numout,*) 'violation hfx  [GW]           (',cd_routine,') = ',zhfx 
    201200      ENDIF 
    202201      ! 
     
    215214      INTEGER  ::   ialert_id         ! number of the current alert 
    216215      REAL(wp) ::   ztmelts           ! ice layer melting point 
    217       CHARACTER (len=30), DIMENSION(20)      ::   cl_alname   ! name of alert 
    218       INTEGER           , DIMENSION(20)      ::   inb_alp     ! number of alerts positive 
     216      CHARACTER (len=30), DIMENSION(20) ::   cl_alname   ! name of alert 
     217      INTEGER           , DIMENSION(20) ::   inb_alp     ! number of alerts positive 
    219218      !!------------------------------------------------------------------- 
    220219 
Note: See TracChangeset for help on using the changeset viewer.