- Timestamp:
- 2017-09-12T20:46:13+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
r8514 r8517 47 47 CONTAINS 48 48 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 ) 50 50 !!---------------------------------------------------------------------- 51 51 !! *** ROUTINE ice_cons_hsm *** … … 56 56 !! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true 57 57 !! It prints in ocean.output if there is a violation of conservation at each time-step 58 !! The thresholds (zv_sill, zs_sill, z h_sill) which determine violations are set to58 !! The thresholds (zv_sill, zs_sill, zt_sill) which determine violations are set to 59 59 !! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 60 60 !! For salt and heat thresholds, ice is considered to have a salinity of 10 … … 63 63 INTEGER , INTENT(in) :: icount ! called at: =0 the begining of the routine, =1 the end 64 64 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,zft68 REAL(wp) ::zvmin, zamin, zamax69 REAL(wp) ::zvtrp, zetrp70 REAL(wp) :: zarea, zv_sill, zs_sill, zh_sill71 REAL(wp), PARAMETER ::zconv = 1.e-9 ! convert W to GW and kg to Mt65 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 72 72 !!---------------------------------------------------------------------- 73 73 ! 74 74 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 ! 75 81 ! ! 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 85 85 ! 86 86 ! ! 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(:,:) ) * zconv90 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 ) * zconv87 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 97 97 98 98 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 99 105 100 106 ! salt flux 101 107 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 102 108 & 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 110 110 111 111 ! heat flux 112 112 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 113 113 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 114 & ) * e1e2t(:,:) ) * zconv - zft_b114 & ) * e1e2t(:,:) ) * zconv - pdiag_ft 115 115 116 116 ! outputs 117 zv i= ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv &118 & - zvi_b ) * r1_rdtice - zfw) * rday119 120 zs mv= ( ( glob_sum( SUM( smv_i * rhoic , dim=3 ) * e1e2t ) * zconv &121 & - zsmv_b) * r1_rdtice + zfs ) * rday122 123 z ei= ( 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 ) & 124 124 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv & 125 & - zei_b) * r1_rdtice + zft125 & - pdiag_t ) * r1_rdtice + zft 126 126 127 127 ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative … … 137 137 zv_sill = zarea * 2.5e-5 138 138 zs_sill = zarea * 25.e-5 139 z h_sill = zarea * 10.e-5139 zt_sill = zarea * 10.e-5 140 140 141 141 IF(lwp) THEN 142 IF ( ABS( zv i ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zvi143 IF ( ABS( zs mv ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zsmv144 IF ( ABS( z ei ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zei142 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 145 145 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 156 153 ENDIF 157 154 ! … … 169 166 !! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true 170 167 !! It prints in ocean.output if there is a violation of conservation at each time-step 171 !! The thresholds (zv_sill, zs_sill, z h_sill) which determine the violation are set to168 !! The thresholds (zv_sill, zs_sill, zt_sill) which determine the violation are set to 172 169 !! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 173 170 !! For salt and heat thresholds, ice is considered to have a salinity of 10 … … 176 173 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 177 174 REAL(wp) :: zhfx, zsfx, zvfx 178 REAL(wp) :: zarea, zv_sill, zs_sill, z h_sill175 REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill 179 176 REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt 180 177 !!---------------------------------------------------------------------- 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 181 184 182 185 ! heat flux … … 184 187 ! & - SUM( qevap_ice * a_i_b, dim=3 ) & !!clem: I think this line must be commented (but need check) 185 188 & ) * e1e2t ) * zconv 186 ! salt flux187 zsfx = glob_sum( ( sfx + diag_smvi ) * e1e2t ) * zconv * rday188 ! water flux189 zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday190 189 191 190 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 193 192 zv_sill = zarea * 2.5e-5 194 193 zs_sill = zarea * 25.e-5 195 z h_sill = zarea * 10.e-5194 zt_sill = zarea * 10.e-5 196 195 197 196 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 ) > z h_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 201 200 ENDIF 202 201 ! … … 215 214 INTEGER :: ialert_id ! number of the current alert 216 215 REAL(wp) :: ztmelts ! ice layer melting point 217 CHARACTER (len=30), DIMENSION(20) 218 INTEGER , DIMENSION(20) 216 CHARACTER (len=30), DIMENSION(20) :: cl_alname ! name of alert 217 INTEGER , DIMENSION(20) :: inb_alp ! number of alerts positive 219 218 !!------------------------------------------------------------------- 220 219
Note: See TracChangeset
for help on using the changeset viewer.
