Changeset 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icectl.F90
- Timestamp:
- 2019-10-29T11:41:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icectl.F90
r10581 r11822 12 12 !! 'key_si3' SI3 sea-ice model 13 13 !!---------------------------------------------------------------------- 14 !! ice_cons_hsm : conservation tests on heat, salt and mass 15 !! ice_cons_final : conservation tests on heat, salt and mass at end of time step 14 !! ice_cons_hsm : conservation tests on heat, salt and mass during a time step (global) 15 !! ice_cons_final : conservation tests on heat, salt and mass at end of time step (global) 16 !! ice_cons2D : conservation tests on heat, salt and mass at each gridcell 16 17 !! ice_ctl : control prints in case of crash 17 18 !! ice_prt : control prints at a given grid point … … 27 28 ! 28 29 USE in_out_manager ! I/O manager 30 USE iom ! I/O manager library 29 31 USE lib_mpp ! MPP library 30 32 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) … … 37 39 PUBLIC ice_cons_hsm 38 40 PUBLIC ice_cons_final 41 PUBLIC ice_cons2D 39 42 PUBLIC ice_ctl 40 43 PUBLIC ice_prt 41 44 PUBLIC ice_prt3D 42 45 46 ! thresold rates for conservation 47 ! these values are changed by the namelist parameter rn_icechk, so that threshold = zchk * rn_icechk 48 REAL(wp), PARAMETER :: zchk_m = 2.5e-7 ! kg/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost 49 REAL(wp), PARAMETER :: zchk_s = 2.5e-6 ! g/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg) 50 REAL(wp), PARAMETER :: zchk_t = 7.5e-2 ! W/m2 <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg) 51 43 52 !! * Substitutions 44 53 # include "vectopt_loop_substitute.h90" … … 59 68 !! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true 60 69 !! It prints in ocean.output if there is a violation of conservation at each time-step 61 !! The thresholds (zv_sill, zs_sill, zt_sill) which determine violations are set to 62 !! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 70 !! The thresholds (zchk_m, zchk_s, zchk_t) determine violations 63 71 !! For salt and heat thresholds, ice is considered to have a salinity of 10 64 72 !! and a heat content of 3e5 J/kg (=latent heat of fusion) … … 68 76 REAL(wp) , INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 69 77 !! 70 REAL(wp) :: z v, zs, zt, zfs, zfv, zft71 REAL(wp) :: zvmin, zamin, zamax78 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, & 79 & zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin 72 80 REAL(wp) :: zvtrp, zetrp 73 REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill 74 REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt 81 REAL(wp) :: zarea 75 82 !!------------------------------------------------------------------- 76 83 ! 77 84 IF( icount == 0 ) THEN 78 ! ! water flux 79 pdiag_fv = glob_sum( 'icectl', & 80 & -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 81 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:) + & 82 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + & 83 & wfx_ice_sub(:,:) + wfx_spr(:,:) & 84 & ) * e1e2t(:,:) ) * zconv 85 86 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) 87 pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) 88 pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) 89 90 ! mass flux 91 pdiag_fv = glob_sum( 'icectl', & 92 & ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 93 & wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) * e1e2t ) 94 ! salt flux 95 pdiag_fs = glob_sum( 'icectl', & 96 & ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + & 97 & sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) 98 ! heat flux 99 pdiag_ft = glob_sum( 'icectl', & 100 & ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 101 & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t ) 102 103 ELSEIF( icount == 1 ) THEN 104 105 ! -- mass diag -- ! 106 zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_rdtice & 107 & + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + & 108 & wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & 109 & wfx_ice_sub + wfx_spr ) * e1e2t ) & 110 & - pdiag_fv 85 111 ! 86 ! ! salt flux87 pdiag_fs = glob_sum( 'icectl',&88 & ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +&89 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)&90 & ) * e1e2t(:,:) ) * zconv112 ! -- salt diag -- ! 113 zdiag_salt = ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) - pdiag_s ) * r1_rdtice & 114 & + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + & 115 & sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) & 116 & - pdiag_fs 91 117 ! 92 ! ! heat flux 93 pdiag_ft = glob_sum( 'icectl', & 94 & ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 95 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 96 & ) * e1e2t(:,:) ) * zconv 97 98 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t * zconv ) 99 100 pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t * zconv ) 101 102 pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 103 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv 104 105 ELSEIF( icount == 1 ) THEN 106 107 ! water flux 108 zfv = glob_sum( 'icectl', & 109 & -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 110 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:) + & 111 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + & 112 & wfx_ice_sub(:,:) + wfx_spr(:,:) & 113 & ) * e1e2t(:,:) ) * zconv - pdiag_fv 114 115 ! salt flux 116 zfs = glob_sum( 'icectl', & 117 & ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 118 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) & 119 & ) * e1e2t(:,:) ) * zconv - pdiag_fs 120 121 ! heat flux 122 zft = glob_sum( 'icectl', & 123 & ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 124 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 125 & ) * e1e2t(:,:) ) * zconv - pdiag_ft 126 127 ! outputs 128 zv = ( ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) * zconv & 129 & - pdiag_v ) * r1_rdtice - zfv ) * rday 130 131 zs = ( ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) * zconv & 132 & - pdiag_s ) * r1_rdtice + zfs ) * rday 133 134 zt = ( glob_sum( 'icectl', & 135 & ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 136 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv & 137 & - pdiag_t ) * r1_rdtice + zft 138 139 ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 140 zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) * zconv * rday 141 zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) * zconv 142 143 zvmin = glob_min( 'icectl', v_i ) 144 zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 145 zamin = glob_min( 'icectl', a_i ) 146 147 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 148 zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2 149 zv_sill = zarea * 2.5e-5 150 zs_sill = zarea * 25.e-5 151 zt_sill = zarea * 10.e-5 152 153 IF(lwp) THEN 154 IF ( ABS( zv ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zv 155 IF ( ABS( zs ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 156 IF ( ABS( zt ) > zt_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zt 157 IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',zvmin 158 IF ( zamax > MAX( rn_amax_n, rn_amax_s ) + epsi10 & 159 & .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' .AND. cd_routine /= 'Hbig' ) & 160 & WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax 161 IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin 162 !clem: the following check fails when using UMx advection scheme (see comments in icedyn_adv.F90) 163 ! IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN 164 ! WRITE(numout,*) 'violation vtrp [Mt/day] (',cd_routine,') = ',zvtrp 165 ! WRITE(numout,*) 'violation etrp [GW] (',cd_routine,') = ',zetrp 166 ! ENDIF 118 ! -- heat diag -- ! 119 zdiag_heat = ( glob_sum( 'icectl', ( SUM(SUM(e_i, dim=4), dim=3) + SUM(SUM(e_s, dim=4), dim=3) ) * e1e2t ) - pdiag_t & 120 & ) * r1_rdtice & 121 & + glob_sum( 'icectl', ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 122 & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t ) & 123 & - pdiag_ft 124 125 ! -- min/max diag -- ! 126 zdiag_amax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 127 zdiag_vmin = glob_min( 'icectl', v_i ) 128 zdiag_amin = glob_min( 'icectl', a_i ) 129 zdiag_smin = glob_min( 'icectl', sv_i ) 130 zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 131 zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 132 133 ! -- advection scheme is conservative? -- ! 134 zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 (only for Prather) 135 zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) ! must be close to 0 (only for Prather) 136 137 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 138 zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 139 140 IF( lwp ) THEN 141 ! check conservation issues 142 IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zarea ) & 143 & WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice 144 IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 145 & WRITE(numout,*) cd_routine,' : violation salt cons. [g] = ',zdiag_salt * rdt_ice 146 IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 147 & WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rdt_ice 148 ! check negative values 149 IF( zdiag_vmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zdiag_vmin 150 IF( zdiag_amin < 0. ) WRITE(numout,*) cd_routine,' : violation a_i < 0 = ',zdiag_amin 151 IF( zdiag_smin < 0. ) WRITE(numout,*) cd_routine,' : violation s_i < 0 = ',zdiag_smin 152 IF( zdiag_eimin < 0. ) WRITE(numout,*) cd_routine,' : violation e_i < 0 = ',zdiag_eimin 153 IF( zdiag_esmin < 0. ) WRITE(numout,*) cd_routine,' : violation e_s < 0 = ',zdiag_esmin 154 ! check maximum ice concentration 155 IF( zdiag_amax > MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 156 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_amax 157 ! check if advection scheme is conservative 158 ! only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc) 159 ! so the formulation for conservation is different (and not coded) 160 ! it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now 161 !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 162 ! & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 167 163 ENDIF 168 164 ! … … 171 167 END SUBROUTINE ice_cons_hsm 172 168 173 174 169 SUBROUTINE ice_cons_final( cd_routine ) 175 170 !!------------------------------------------------------------------- … … 180 175 !! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true 181 176 !! It prints in ocean.output if there is a violation of conservation at each time-step 182 !! The thresholds (zv_sill, zs_sill, zt_sill) which determine the violation are set to 183 !! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 177 !! The thresholds (zchk_m, zchk_s, zchk_t) determine the violations 184 178 !! For salt and heat thresholds, ice is considered to have a salinity of 10 185 179 !! and a heat content of 3e5 J/kg (=latent heat of fusion) 186 180 !!------------------------------------------------------------------- 187 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 188 REAL(wp) :: zqmass, zhfx, zsfx, zvfx 189 REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill 190 REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt 181 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 182 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat 183 REAL(wp) :: zarea 191 184 !!------------------------------------------------------------------- 192 185 193 186 ! water flux 194 zvfx = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday 195 196 ! salt flux 197 zsfx = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) * zconv * rday 198 199 ! heat flux 187 ! -- mass diag -- ! 188 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 189 190 ! -- salt diag -- ! 191 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) 192 193 ! -- heat diag -- ! 200 194 ! clem: not the good formulation 201 !!zhfx = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr & 202 !! & ) * e1e2t ) * zconv 203 204 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 205 zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2 206 zv_sill = zarea * 2.5e-5 207 zs_sill = zarea * 25.e-5 208 zt_sill = zarea * 10.e-5 209 210 IF(lwp) THEN 211 IF( ABS( zvfx ) > zv_sill ) WRITE(numout,*) 'violation vfx [Mt/day] (',cd_routine,') = ',zvfx 212 IF( ABS( zsfx ) > zs_sill ) WRITE(numout,*) 'violation sfx [psu*Mt/day] (',cd_routine,') = ',zsfx 213 !!IF( ABS( zhfx ) > zt_sill ) WRITE(numout,*) 'violation hfx [GW] (',cd_routine,') = ',zhfx 195 !!zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr & 196 !! & ) * e1e2t ) 197 198 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 199 zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 200 201 IF( lwp ) THEN 202 IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zarea ) & 203 & WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice 204 IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 205 & WRITE(numout,*) cd_routine,' : violation salt cons. [g] = ',zdiag_salt * rdt_ice 206 !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rdt_ice 214 207 ENDIF 215 208 ! 216 209 END SUBROUTINE ice_cons_final 217 210 211 SUBROUTINE ice_cons2D( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft ) 212 !!------------------------------------------------------------------- 213 !! *** ROUTINE ice_cons2D *** 214 !! 215 !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine 216 !! + test if ice concentration and volume are > 0 217 !! 218 !! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true 219 !! It stops the code if there is a violation of conservation at any gridcell 220 !!------------------------------------------------------------------- 221 INTEGER , INTENT(in) :: icount ! called at: =0 the begining of the routine, =1 the end 222 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 223 REAL(wp) , DIMENSION(jpi,jpj), INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 224 !! 225 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass, zdiag_salt, zdiag_heat, & 226 & zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax 227 INTEGER :: jl, jk 228 LOGICAL :: ll_stop_m = .FALSE. 229 LOGICAL :: ll_stop_s = .FALSE. 230 LOGICAL :: ll_stop_t = .FALSE. 231 CHARACTER(len=120) :: clnam ! filename for the output 232 !!------------------------------------------------------------------- 233 ! 234 IF( icount == 0 ) THEN 235 236 pdiag_v = SUM( v_i * rhoi + v_s * rhos, dim=3 ) 237 pdiag_s = SUM( sv_i * rhoi , dim=3 ) 238 pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 239 240 ! mass flux 241 pdiag_fv = wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 242 & wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr 243 ! salt flux 244 pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam 245 ! heat flux 246 pdiag_ft = hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 247 & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr 248 249 ELSEIF( icount == 1 ) THEN 250 251 ! -- mass diag -- ! 252 zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos, dim=3 ) - pdiag_v ) * r1_rdtice & 253 & + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 254 & wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) & 255 & - pdiag_fv 256 IF( MAXVAL( ABS(zdiag_mass) ) > zchk_m * rn_icechk_cel ) ll_stop_m = .TRUE. 257 ! 258 ! -- salt diag -- ! 259 zdiag_salt = ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_rdtice & 260 & + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) & 261 & - pdiag_fs 262 IF( MAXVAL( ABS(zdiag_salt) ) > zchk_s * rn_icechk_cel ) ll_stop_s = .TRUE. 263 ! 264 ! -- heat diag -- ! 265 zdiag_heat = ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) - pdiag_t ) * r1_rdtice & 266 & + ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 267 & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) & 268 & - pdiag_ft 269 IF( MAXVAL( ABS(zdiag_heat) ) > zchk_t * rn_icechk_cel ) ll_stop_t = .TRUE. 270 ! 271 ! -- other diags -- ! 272 ! a_i < 0 273 zdiag_amin(:,:) = 0._wp 274 DO jl = 1, jpl 275 WHERE( a_i(:,:,jl) < 0._wp ) zdiag_amin(:,:) = 1._wp 276 ENDDO 277 ! v_i < 0 278 zdiag_vmin(:,:) = 0._wp 279 DO jl = 1, jpl 280 WHERE( v_i(:,:,jl) < 0._wp ) zdiag_vmin(:,:) = 1._wp 281 ENDDO 282 ! s_i < 0 283 zdiag_smin(:,:) = 0._wp 284 DO jl = 1, jpl 285 WHERE( s_i(:,:,jl) < 0._wp ) zdiag_smin(:,:) = 1._wp 286 ENDDO 287 ! e_i < 0 288 zdiag_emin(:,:) = 0._wp 289 DO jl = 1, jpl 290 DO jk = 1, nlay_i 291 WHERE( e_i(:,:,jk,jl) < 0._wp ) zdiag_emin(:,:) = 1._wp 292 ENDDO 293 ENDDO 294 ! a_i > amax 295 !WHERE( SUM( a_i, dim=3 ) > ( MAX( rn_amax_n, rn_amax_s ) + epsi10 ) ; zdiag_amax(:,:) = 1._wp 296 !ELSEWHERE ; zdiag_amax(:,:) = 0._wp 297 !END WHERE 298 299 IF( ll_stop_m .OR. ll_stop_s .OR. ll_stop_t ) THEN 300 clnam = 'diag_ice_conservation_'//cd_routine 301 CALL ice_cons_wri( clnam, zdiag_mass, zdiag_salt, zdiag_heat, zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin ) 302 ENDIF 303 304 IF( ll_stop_m ) CALL ctl_stop( 'STOP', cd_routine//': ice mass conservation issue' ) 305 IF( ll_stop_s ) CALL ctl_stop( 'STOP', cd_routine//': ice salt conservation issue' ) 306 IF( ll_stop_t ) CALL ctl_stop( 'STOP', cd_routine//': ice heat conservation issue' ) 307 308 ENDIF 309 310 END SUBROUTINE ice_cons2D 311 312 SUBROUTINE ice_cons_wri( cdfile_name, pdiag_mass, pdiag_salt, pdiag_heat, pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin ) 313 !!--------------------------------------------------------------------- 314 !! *** ROUTINE ice_cons_wri *** 315 !! 316 !! ** Purpose : create a NetCDF file named cdfile_name which contains 317 !! the instantaneous fields when conservation issue occurs 318 !! 319 !! ** Method : NetCDF files using ioipsl 320 !!---------------------------------------------------------------------- 321 CHARACTER(len=*), INTENT( in ) :: cdfile_name ! name of the file created 322 REAL(wp), DIMENSION(:,:), INTENT( in ) :: pdiag_mass, pdiag_salt, pdiag_heat, & 323 & pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin !!, pdiag_amax 324 !! 325 INTEGER :: inum 326 !!---------------------------------------------------------------------- 327 ! 328 IF(lwp) WRITE(numout,*) 329 IF(lwp) WRITE(numout,*) 'ice_cons_wri : single instantaneous ice state' 330 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ named :', cdfile_name, '...nc' 331 IF(lwp) WRITE(numout,*) 332 333 CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 334 335 CALL iom_rstput( 0, 0, inum, 'cons_mass', pdiag_mass(:,:) , ktype = jp_r8 ) ! ice mass spurious lost/gain 336 CALL iom_rstput( 0, 0, inum, 'cons_salt', pdiag_salt(:,:) , ktype = jp_r8 ) ! ice salt spurious lost/gain 337 CALL iom_rstput( 0, 0, inum, 'cons_heat', pdiag_heat(:,:) , ktype = jp_r8 ) ! ice heat spurious lost/gain 338 ! other diags 339 CALL iom_rstput( 0, 0, inum, 'aneg_count', pdiag_amin(:,:) , ktype = jp_r8 ) ! 340 CALL iom_rstput( 0, 0, inum, 'vneg_count', pdiag_vmin(:,:) , ktype = jp_r8 ) ! 341 CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 ) ! 342 CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 ) ! 343 344 CALL iom_close( inum ) 345 346 END SUBROUTINE ice_cons_wri 218 347 219 348 SUBROUTINE ice_ctl( kt ) … … 238 367 ialert_id = 2 ! reference number of this alert 239 368 cl_alname(ialert_id) = ' Incompat vol and con ' ! name of the alert 240 241 369 DO jl = 1, jpl 242 370 DO jj = 1, jpj 243 371 DO ji = 1, jpi 244 372 IF( v_i(ji,jj,jl) /= 0._wp .AND. a_i(ji,jj,jl) == 0._wp ) THEN 245 !WRITE(numout,*) ' ALERTE 2 : Incompatible volume and concentration ' 246 !WRITE(numout,*) ' at_i ', at_i(ji,jj) 247 !WRITE(numout,*) ' Point - category', ji, jj, jl 248 !WRITE(numout,*) ' a_i *** a_i_b ', a_i (ji,jj,jl), a_i_b (ji,jj,jl) 249 !WRITE(numout,*) ' v_i *** v_i_b ', v_i (ji,jj,jl), v_i_b (ji,jj,jl) 373 WRITE(numout,*) ' ALERTE 2 : Incompatible volume and concentration ' 250 374 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 251 375 ENDIF … … 261 385 DO ji = 1, jpi 262 386 IF( h_i(ji,jj,jl) > 50._wp ) THEN 387 WRITE(numout,*) ' ALERTE 3 : Very thick ice' 263 388 !CALL ice_prt( kt, ji, jj, 2, ' ALERTE 3 : Very thick ice ' ) 264 389 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 … … 272 397 DO jj = 1, jpj 273 398 DO ji = 1, jpi 274 IF( MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5.AND. &399 IF( MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2. .AND. & 275 400 & at_i(ji,jj) > 0._wp ) THEN 401 WRITE(numout,*) ' ALERTE 4 : Very fast ice' 276 402 !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 4 : Very fast ice ' ) 277 !WRITE(numout,*) ' ice strength : ', strength(ji,jj) 278 !WRITE(numout,*) ' oceanic stress utau : ', utau(ji,jj) 279 !WRITE(numout,*) ' oceanic stress vtau : ', vtau(ji,jj) 280 !WRITE(numout,*) ' sea-ice stress utau_ice : ', utau_ice(ji,jj) 281 !WRITE(numout,*) ' sea-ice stress vtau_ice : ', vtau_ice(ji,jj) 282 !WRITE(numout,*) ' sst : ', sst_m(ji,jj) 283 !WRITE(numout,*) ' sss : ', sss_m(ji,jj) 284 !WRITE(numout,*) 403 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 404 ENDIF 405 END DO 406 END DO 407 408 ! Alert on salt flux 409 ialert_id = 5 ! reference number of this alert 410 cl_alname(ialert_id) = ' High salt flux ' ! name of the alert 411 DO jj = 1, jpj 412 DO ji = 1, jpi 413 IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN ! = 1 psu/day for 1m ocean depth 414 WRITE(numout,*) ' ALERTE 5 : High salt flux' 415 !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 : High salt flux ' ) 285 416 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 286 417 ENDIF … … 294 425 DO ji = 1, jpi 295 426 IF( tmask(ji,jj,1) <= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 427 WRITE(numout,*) ' ALERTE 6 : Ice on continents' 296 428 !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 6 : Ice on continents ' ) 297 !WRITE(numout,*) ' masks s, u, v : ', tmask(ji,jj,1), umask(ji,jj,1), vmask(ji,jj,1)298 !WRITE(numout,*) ' sst : ', sst_m(ji,jj)299 !WRITE(numout,*) ' sss : ', sss_m(ji,jj)300 !WRITE(numout,*) ' at_i(ji,jj) : ', at_i(ji,jj)301 !WRITE(numout,*) ' v_ice(ji,jj) : ', v_ice(ji,jj)302 !WRITE(numout,*) ' v_ice(ji,jj-1) : ', v_ice(ji,jj-1)303 !WRITE(numout,*) ' u_ice(ji-1,jj) : ', u_ice(ji-1,jj)304 !WRITE(numout,*) ' u_ice(ji,jj) : ', v_ice(ji,jj)305 !306 429 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 307 430 ENDIF … … 317 440 DO ji = 1, jpi 318 441 IF( s_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 442 WRITE(numout,*) ' ALERTE 7 : Very fresh ice' 319 443 ! CALL ice_prt(kt,ji,jj,1, ' ALERTE 7 : Very fresh ice ' ) 320 ! WRITE(numout,*) ' sst : ', sst_m(ji,jj)321 ! WRITE(numout,*) ' sss : ', sss_m(ji,jj)322 ! WRITE(numout,*)323 444 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 324 445 ENDIF … … 327 448 END DO 328 449 ! 450 ! Alert if qns very big 451 ialert_id = 8 ! reference number of this alert 452 cl_alname(ialert_id) = ' fnsolar very big ' ! name of the alert 453 DO jj = 1, jpj 454 DO ji = 1, jpi 455 IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 456 ! 457 WRITE(numout,*) ' ALERTE 8 : Very high non-solar heat flux' 458 !CALL ice_prt( kt, ji, jj, 2, ' ') 459 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 460 ! 461 ENDIF 462 END DO 463 END DO 464 !+++++ 329 465 330 466 ! ! Alert if too old ice … … 337 473 ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 338 474 ( a_i(ji,jj,jl) > 0._wp ) ) THEN 475 WRITE(numout,*) ' ALERTE 9 : Wrong ice age' 339 476 !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 9 : Wrong ice age ') 340 477 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 … … 343 480 END DO 344 481 END DO 345 346 ! Alert on salt flux 347 ialert_id = 5 ! reference number of this alert 348 cl_alname(ialert_id) = ' High salt flux ' ! name of the alert 349 DO jj = 1, jpj 350 DO ji = 1, jpi 351 IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN ! = 1 psu/day for 1m ocean depth 352 !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 : High salt flux ' ) 353 !DO jl = 1, jpl 354 !WRITE(numout,*) ' Category no: ', jl 355 !WRITE(numout,*) ' a_i : ', a_i (ji,jj,jl) , ' a_i_b : ', a_i_b (ji,jj,jl) 356 !WRITE(numout,*) ' v_i : ', v_i (ji,jj,jl) , ' v_i_b : ', v_i_b (ji,jj,jl) 357 !WRITE(numout,*) ' ' 358 !END DO 359 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 360 ENDIF 361 END DO 362 END DO 363 364 ! Alert if qns very big 365 ialert_id = 8 ! reference number of this alert 366 cl_alname(ialert_id) = ' fnsolar very big ' ! name of the alert 367 DO jj = 1, jpj 368 DO ji = 1, jpi 369 IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 370 ! 371 !WRITE(numout,*) ' ALERTE 8 : Very high non-solar heat flux' 372 !WRITE(numout,*) ' ji, jj : ', ji, jj 373 !WRITE(numout,*) ' qns : ', qns(ji,jj) 374 !WRITE(numout,*) ' sst : ', sst_m(ji,jj) 375 !WRITE(numout,*) ' sss : ', sss_m(ji,jj) 376 ! 377 !CALL ice_prt( kt, ji, jj, 2, ' ') 378 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 379 ! 380 ENDIF 381 END DO 382 END DO 383 !+++++ 384 482 385 483 ! Alert if very warm ice 386 484 ialert_id = 10 ! reference number of this alert … … 392 490 DO ji = 1, jpi 393 491 ztmelts = -rTmlt * sz_i(ji,jj,jk,jl) + rt0 394 IF( t_i(ji,jj,jk,jl) >= ztmelts .AND. v_i(ji,jj,jl) > 1.e-10 & 395 & .AND. a_i(ji,jj,jl) > 0._wp ) THEN 396 !WRITE(numout,*) ' ALERTE 10 : Very warm ice' 397 !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 398 !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl) 399 !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl) 400 !WRITE(numout,*) ' sz_i: ', sz_i(ji,jj,jk,jl) 401 !WRITE(numout,*) ' ztmelts : ', ztmelts 402 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 492 IF( t_i(ji,jj,jk,jl) > ztmelts .AND. v_i(ji,jj,jl) > 1.e-10 & 493 & .AND. a_i(ji,jj,jl) > 0._wp ) THEN 494 WRITE(numout,*) ' ALERTE 10 : Very warm ice' 495 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 403 496 ENDIF 404 497 END DO … … 427 520 END SUBROUTINE ice_ctl 428 521 429 430 522 SUBROUTINE ice_prt( kt, ki, kj, kn, cd1 ) 431 523 !!------------------------------------------------------------------- … … 435 527 !! in ocean.ouput 436 528 !! 3 possibilities exist 437 !! n = 1/-1 -> simple ice state (plus Mechanical Check if -1)529 !! n = 1/-1 -> simple ice state 438 530 !! n = 2 -> exhaustive state 439 531 !! n = 3 -> ice/ocean salt fluxes … … 474 566 WRITE(numout,*) ' - Cell values ' 475 567 WRITE(numout,*) ' ~~~~~~~~~~~ ' 476 WRITE(numout,*) ' cell area : ', e1e2t(ji,jj)477 568 WRITE(numout,*) ' at_i : ', at_i(ji,jj) 569 WRITE(numout,*) ' ato_i : ', ato_i(ji,jj) 478 570 WRITE(numout,*) ' vt_i : ', vt_i(ji,jj) 479 571 WRITE(numout,*) ' vt_s : ', vt_s(ji,jj) … … 495 587 END DO 496 588 ENDIF 497 IF( kn == -1 ) THEN498 WRITE(numout,*) ' Mechanical Check ************** '499 WRITE(numout,*) ' Check what means ice divergence '500 WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj)501 WRITE(numout,*) ' Total lead fraction ', ato_i(ji,jj)502 WRITE(numout,*) ' Sum of both ', ato_i(ji,jj) + at_i(ji,jj)503 WRITE(numout,*) ' Sum of both minus 1 ', ato_i(ji,jj) + at_i(ji,jj) - 1.00504 ENDIF505 506 589 507 590 !-------------------- … … 517 600 WRITE(numout,*) ' - Cell values ' 518 601 WRITE(numout,*) ' ~~~~~~~~~~~ ' 519 WRITE(numout,*) ' cell area : ', e1e2t(ji,jj)520 602 WRITE(numout,*) ' at_i : ', at_i(ji,jj) 521 603 WRITE(numout,*) ' vt_i : ', vt_i(ji,jj) … … 616 698 !! 617 699 !!------------------------------------------------------------------- 618 CHARACTER(len=*), INTENT(in) ::cd_routine ! name of the routine619 INTEGER ::jk, jl ! dummy loop indices700 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 701 INTEGER :: jk, jl ! dummy loop indices 620 702 621 703 CALL prt_ctl_info(' ========== ') … … 676 758 677 759 END SUBROUTINE ice_prt3D 678 760 679 761 #else 680 762 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.