- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icethd_do.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/ICE/icethd_do.F90
r12178 r12928 44 44 REAL(wp) :: rn_Cfraz ! squeezing coefficient for collection of bottom frazil ice 45 45 46 !! * Substitutions 47 # include "do_loop_substitute.h90" 46 48 !!---------------------------------------------------------------------- 47 49 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 139 141 ! Physical constants 140 142 zhicrit = 0.04 ! frazil ice thickness 141 ztwogp = 2. * r au0 / ( grav * 0.3 * ( rau0 - rhoi ) ) ! reduced grav143 ztwogp = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoi ) ) ! reduced grav 142 144 zsqcd = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag) 143 145 zgamafr = 0.03 144 146 ! 145 DO jj = 2, jpjm1 146 DO ji = 2, jpim1 147 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast 148 ! -- Wind stress -- ! 149 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & 150 & + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp 151 ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) & 152 & + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp 153 ! Square root of wind stress 154 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 155 156 ! -- Frazil ice velocity -- ! 157 rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 158 zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 159 zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 160 161 ! -- Pack ice velocity -- ! 162 zvgx = ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 163 zvgy = ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 164 165 ! -- Relative frazil/pack ice velocity -- ! 166 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 167 zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) & 168 & + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch 169 zvrel(ji,jj) = SQRT( zvrel2 ) 170 171 ! -- new ice thickness (iterative loop) -- ! 172 ht_i_new(ji,jj) = zhicrit + ( zhicrit + 0.1 ) & 173 & / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) - zhicrit * zhicrit ) * ztwogp * zvrel2 174 175 iter = 1 176 DO WHILE ( iter < 20 ) 177 zf = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) - & 178 & ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2 179 zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0 * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 180 181 ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 ) 182 iter = iter + 1 183 END DO 184 ! 185 ! bound ht_i_new (though I don't see why it should be necessary) 186 ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) ) 187 ! 188 ENDIF 147 DO_2D_00_00 148 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast 149 ! -- Wind stress -- ! 150 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & 151 & + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp 152 ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) & 153 & + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp 154 ! Square root of wind stress 155 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 156 157 ! -- Frazil ice velocity -- ! 158 rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 159 zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 160 zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 161 162 ! -- Pack ice velocity -- ! 163 zvgx = ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 164 zvgy = ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 165 166 ! -- Relative frazil/pack ice velocity -- ! 167 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 168 zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) & 169 & + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch 170 zvrel(ji,jj) = SQRT( zvrel2 ) 171 172 ! -- new ice thickness (iterative loop) -- ! 173 ht_i_new(ji,jj) = zhicrit + ( zhicrit + 0.1 ) & 174 & / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) - zhicrit * zhicrit ) * ztwogp * zvrel2 175 176 iter = 1 177 DO WHILE ( iter < 20 ) 178 zf = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) - & 179 & ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2 180 zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0 * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 181 182 ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 ) 183 iter = iter + 1 184 END DO 189 185 ! 190 END DO 191 END DO 186 ! bound ht_i_new (though I don't see why it should be necessary) 187 ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) ) 188 ! 189 ENDIF 190 ! 191 END_2D 192 192 ! 193 193 CALL lbc_lnk_multi( 'icethd_do', zvrel, 'T', 1., ht_i_new, 'T', 1. ) … … 202 202 ! Identify grid points where new ice forms 203 203 npti = 0 ; nptidx(:) = 0 204 DO jj = 1, jpj 205 DO ji = 1, jpi 206 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN 207 npti = npti + 1 208 nptidx( npti ) = (jj - 1) * jpi + ji 209 ENDIF 210 END DO 211 END DO 204 DO_2D_11_11 205 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN 206 npti = npti + 1 207 nptidx( npti ) = (jj - 1) * jpi + ji 208 ENDIF 209 END_2D 212 210 213 211 ! Move from 2-D to 1-D vectors … … 291 289 292 290 ! Contribution to heat flux to the ocean [W.m-2], >0 293 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_ rdtice291 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_Dt_ice 294 292 ! Total heat flux used in this process [W.m-2] 295 hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_ rdtice293 hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_Dt_ice 296 294 ! mass flux 297 wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_ rdtice295 wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_Dt_ice 298 296 ! salt flux 299 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_ rdtice297 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_Dt_ice 300 298 END DO 301 299 … … 443 441 !!------------------------------------------------------------------- 444 442 ! 445 REWIND( numnam_ice_ref ) ! Namelist namthd_do in reference namelist : Ice thermodynamics446 443 READ ( numnam_ice_ref, namthd_do, IOSTAT = ios, ERR = 901) 447 444 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_do in reference namelist' ) 448 REWIND( numnam_ice_cfg ) ! Namelist namthd_do in configuration namelist : Ice thermodynamics449 445 READ ( numnam_ice_cfg, namthd_do, IOSTAT = ios, ERR = 902 ) 450 446 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namthd_do in configuration namelist' )
Note: See TracChangeset
for help on using the changeset viewer.
