Changeset 5070
- Timestamp:
- 2015-02-09T14:39:07+01:00 (9 years ago)
- Location:
- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO
- Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r5067 r5070 147 147 !! tm_i | - | Mean sea ice temperature | K | 148 148 !! ot_i ! - ! Sea ice areal age content | day | 149 !! et_i ! - ! Total ice enthalpy | 10^9 J|150 !! et_s ! - ! Total snow enthalpy | 10^9 J|149 !! et_i ! - ! Total ice enthalpy | J/m2 | 150 !! et_s ! - ! Total snow enthalpy | J/m2 | 151 151 !! bv_i ! - ! Mean relative brine volume | ??? | 152 152 !!===================================================================== … … 172 172 173 173 ! !!** ice-dynamics namelist (namicedyn) ** 174 LOGICAL , PUBLIC :: ln_icestr_bvf !: use brine volume to diminish ice strength 174 175 INTEGER , PUBLIC :: nn_icestr !: ice strength parameterization (0=Hibler79 1=Rothrock75) 175 INTEGER , PUBLIC :: nn_icestr_bvf !: use brine volume to diminish ice strength176 176 INTEGER , PUBLIC :: nn_nevp !: number of iterations for subcycling 177 INTEGER , PUBLIC :: nn_ahi0 !: sea-ice hor. eddy diffusivity coeff. (3 ways of calculation) 178 REAL(wp), PUBLIC :: rn_pe_rdg !: ridging work divided by pot. energy change in ridging, nn_icestr = 1 177 179 REAL(wp), PUBLIC :: rn_cio !: drag coefficient for oceanic stress 178 180 REAL(wp), PUBLIC :: rn_pstar !: determines ice strength (N/M), Hibler JPO79 … … 202 204 ! !!** ice-mechanical redistribution namelist (namiceitdme) 203 205 REAL(wp), PUBLIC :: rn_cs !: fraction of shearing energy contributing to ridging 204 REAL(wp), PUBLIC :: rn_pe_rdg !: ridging work divided by pot. energy change in ridging, nn_ice str = 1205 206 REAL(wp), PUBLIC :: rn_fsnowrdg !: fractional snow loss to the ocean during ridging 206 207 REAL(wp), PUBLIC :: rn_fsnowrft !: fractional snow loss to the ocean during ridging … … 217 218 218 219 ! !!** ice-mechanical redistribution namelist (namiceitdme) 219 INTEGER , PUBLIC :: nn_rafting !: rafting of ice or not220 LOGICAL , PUBLIC :: ln_rafting !: rafting of ice or not 220 221 INTEGER , PUBLIC :: nn_partfun !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007) 221 222 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r5067 r5070 173 173 IF( icount == 0 ) THEN 174 174 175 zvi_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 175 zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 176 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 177 & ) * e12t(:,:) * tmask(:,:,1) ) 178 179 zfw_b = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 180 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 181 & ) * e12t(:,:) * tmask(:,:,1) ) 182 183 zft_b = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 184 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 185 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) 186 187 zvi_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 188 176 189 zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 177 zei_b = glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv + & 178 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv ) 179 zfw_b = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 180 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 181 & ) * e12t(:,:) * tmask(:,:,1) ) 182 zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 183 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 184 & ) * e12t(:,:) * tmask(:,:,1) ) 185 zft_b = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 186 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 187 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) 190 191 zei_b = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + & 192 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) & 193 ) * e12t(:,:) * tmask(:,:,1) * zconv ) 188 194 189 195 ELSEIF( icount == 1 ) THEN 190 196 191 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 192 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 193 & ) * e12t(:,:) * tmask(:,:,1) ) - zfs_b 194 zfw = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 195 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 196 & ) * e12t(:,:) * tmask(:,:,1) ) - zfw_b 197 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 198 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 199 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b 197 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 198 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 199 & ) * e12t(:,:) * tmask(:,:,1) ) - zfs_b 200 201 zfw = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 202 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 203 & ) * e12t(:,:) * tmask(:,:,1) ) - zfw_b 204 205 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 206 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 207 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b 200 208 201 zvi = ( glob_sum( SUM( 209 zvi = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) & 202 210 & * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw 211 203 212 zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 204 zei = glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv + & 205 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv ) & 206 & * r1_rdtice - zei_b * r1_rdtice + zft 207 208 zvmin = glob_min(v_i) 209 zamax = glob_max(SUM(a_i,dim=3)) 210 zamin = glob_min(a_i) 213 214 zei = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + & 215 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) & 216 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft 217 218 zvmin = glob_min( v_i ) 219 zamax = glob_max( SUM( a_i, dim=3 ) ) 220 zamin = glob_min( a_i ) 211 221 212 222 IF(lwp) THEN 213 223 IF ( ABS( zvi ) > 1.e-4 ) WRITE(numout,*) 'violation volume [kg/day] (',cd_routine,') = ',(zvi * rday) 214 224 IF ( ABS( zsmv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday) 215 IF ( ABS( zei ) > 1 .e-2) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',(zei)225 IF ( ABS( zei ) > 1 ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',(zei) 216 226 IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',(zvmin) 217 227 IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > rn_amax+epsi10 ) THEN -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
r5067 r5070 180 180 DO jj = 1, jpj 181 181 DO ji = 1, jpi 182 IF( ABS( sfx (ji,jj) ) .GT.1.0e-2 ) THEN ! = 1 psu/day for 1m ocean depth182 IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN ! = 1 psu/day for 1m ocean depth 183 183 !CALL lim_prt( kt, ji, jj, 3, ' ALERTE 5 : High salt flux ' ) 184 184 !DO jl = 1, jpl -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
r5067 r5070 113 113 114 114 ! Heat budget 115 zbg_ihc = glob_sum( et_i(:,:) * 1.e-20 ) ! ice heat content [1.e-20 J]116 zbg_shc = glob_sum( et_s(:,:) * 1.e-20 ) ! snow heat content [1.e-20 J]115 zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content [1.e20 J] 116 zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J] 117 117 zbg_hfx_dhc = glob_sum( diag_heat_dhc(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 118 118 zbg_hfx_spr = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] … … 137 137 z_frc_sal = r1_rau0 * glob_sum( sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) ! salt fluxes 138 138 z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 139 & wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 139 & wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + & 140 & wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 140 141 ! 141 142 frc_vol = frc_vol + z_frc_vol * rdt_ice -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r5067 r5070 172 172 ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean 173 173 ! 174 zcoef = SQRT( 0.5_wp ) /rau0174 zcoef = SQRT( 0.5_wp ) * r1_rau0 175 175 DO jj = 2, jpjm1 176 176 DO ji = fs_2, fs_jpim1 ! vector opt. … … 243 243 !!------------------------------------------------------------------- 244 244 INTEGER :: ios ! Local integer output status for namelist read 245 NAMELIST/namicedyn/ nn_icestr, nn_icestr_bvf, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, nn_nevp, rn_relast, rn_ahi0_ref 245 NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, & 246 & nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref 246 247 INTEGER :: ji, jj 247 248 REAL(wp) :: za00, zd_max … … 261 262 WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 262 263 WRITE(numout,*) '~~~~~~~~~~~~' 263 WRITE(numout,*)' ice strength parameterization (0=Hibler 1=Rothrock) nn_icestr =', nn_icestr 264 WRITE(numout,*)' Including brine volume in ice strength comp. nn_icestr_bvf=', nn_icestr_bvf 265 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 266 WRITE(numout,*) ' first bulk-rheology parameter rn_pstar = ', rn_pstar 267 WRITE(numout,*) ' second bulk-rhelogy parameter rn_crhg = ', rn_crhg 268 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 269 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 270 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 271 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 272 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref 264 WRITE(numout,*)' ice strength parameterization (0=Hibler 1=Rothrock) nn_icestr = ', nn_icestr 265 WRITE(numout,*)' Including brine volume in ice strength comp. ln_icestr_bvf = ', ln_icestr_bvf 266 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging rn_pe_rdg = ', rn_pe_rdg 267 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 268 WRITE(numout,*) ' first bulk-rheology parameter rn_pstar = ', rn_pstar 269 WRITE(numout,*) ' second bulk-rhelogy parameter rn_crhg = ', rn_crhg 270 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 271 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 272 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 273 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 274 WRITE(numout,*) ' horizontal diffusivity calculation nn_ahi0 = ', nn_ahi0 275 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref 273 276 ENDIF 274 277 ! … … 277 280 ! 278 281 ! Diffusion coefficients 279 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 280 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 281 282 za00 = rn_ahi0_ref / ( 1.e05_wp ) ! 1.e05 = 100km = max grid space at 60° latitude in orca2 283 ! (60° = min latitude for ice cover) 284 DO jj = 1, jpj 285 DO ji = 1, jpi 286 ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 287 ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 282 SELECT CASE( nn_ahi0 ) 283 284 CASE( 0 ) 285 ahiu(:,:) = rn_ahi0_ref 286 ahiv(:,:) = rn_ahi0_ref 287 288 IF(lwp) WRITE(numout,*) '' 289 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim constant = rn_ahi0_ref' 290 291 CASE( 1 ) 292 293 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 294 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 295 296 ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2 297 ! (60° = min latitude for ice cover) 298 ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 299 300 IF(lwp) WRITE(numout,*) '' 301 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 302 IF(lwp) WRITE(numout,*) ' value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 303 304 CASE( 2 ) 305 306 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 307 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 308 309 za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2 310 ! (60° = min latitude for ice cover) 311 DO jj = 1, jpj 312 DO ji = 1, jpi 313 ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 314 ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 315 END DO 288 316 END DO 289 END DO 290 ! 291 IF(lwp) WRITE(numout,*) '' 292 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1' 293 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 294 317 ! 318 IF(lwp) WRITE(numout,*) '' 319 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1' 320 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 321 322 END SELECT 323 295 324 END SUBROUTINE lim_dyn_init 296 325 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r5047 r5070 72 72 DO jj = 2, jpjm1 73 73 DO ji = fs_2 , fs_jpim1 ! vector opt. 74 efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj))74 efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj) 75 75 END DO 76 76 END DO … … 98 98 DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction 99 99 DO ji = 1 , fs_jpim1 ! vector opt. 100 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) /e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )101 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) /e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )100 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 101 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 102 102 END DO 103 103 END DO … … 105 105 DO jj= 2, jpjm1 ! diffusive trend : divergence of the fluxes 106 106 DO ji = fs_2 , fs_jpim1 ! vector opt. 107 zdiv (ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj ) & 108 & + zflv(ji,jj) - zflv(ji ,jj-1) ) / ( e1t (ji,jj) * e2t (ji,jj) ) 107 zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj) 109 108 END DO 110 109 END DO … … 116 115 zrlxint = ( ztab0(ji,jj) & 117 116 & + rdt_ice * ( zalfa * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) ) & 118 & + ( 1.0 - zalfa ) * zdiv0(ji,jj) ) )&119 & / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )117 & + ( 1.0 - zalfa ) * zdiv0(ji,jj) ) & 118 & ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) ) 120 119 zrlx(ji,jj) = ptab(ji,jj) + zrelax * ( zrlxint - ptab(ji,jj) ) 121 120 END DO … … 139 138 DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction 140 139 DO ji = 1 , fs_jpim1 ! vector opt. 141 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) /e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )142 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) /e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )140 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 141 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 143 142 END DO 144 143 END DO … … 146 145 DO jj= 2, jpjm1 ! diffusive trend : divergence of the fluxes 147 146 DO ji = fs_2 , fs_jpim1 ! vector opt. 148 zdiv (ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj ) & 149 & + zflv(ji,jj) - zflv(ji ,jj-1) ) / ( e1t (ji,jj) * e2t (ji,jj) ) 147 zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj) 150 148 ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) ) 151 149 END DO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r5067 r5070 157 157 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 158 158 !-----------------------------------------------------------------------------! 159 Cp = 0.5 * grav * (rau0-rhoic) * rhoic /rau0 ! proport const for PE159 Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE 160 160 ! 161 161 CALL lim_itd_me_ridgeprep ! prepare ridging … … 191 191 ! (thick, newly ridged ice). 192 192 193 closing_net(ji,jj) = rn_cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )193 closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 194 194 195 195 ! 2.2 divu_adv … … 235 235 ! Reduce the closing rate if more than 100% of the open water 236 236 ! would be removed. Reduce the opening rate proportionately. 237 IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT.0.0 ) THEN237 IF ( ato_i(ji,jj) > epsi10 .AND. athorn(ji,jj,0) > 0.0 ) THEN 238 238 w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 239 IF ( w1 .GT.ato_i(ji,jj)) THEN239 IF ( w1 > ato_i(ji,jj)) THEN 240 240 tmpfac = ato_i(ji,jj) / w1 241 241 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac … … 286 286 DO jj = 1, jpj 287 287 DO ji = 1, jpi 288 IF (ABS(asum(ji,jj) - kamax ) .LT.epsi10) THEN288 IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN 289 289 closing_net(ji,jj) = 0._wp 290 290 opning (ji,jj) = 0._wp … … 357 357 WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 358 358 END DO 359 ENDIF ! asum 360 361 END DO !ji 362 END DO !jj 359 ENDIF 360 END DO 361 END DO 363 362 364 363 ! Conservation check … … 482 481 * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) ) 483 482 !!gm Optimization: (a**3-b**3)/(a-b) = a*a+ab+b*b ==> less costly operations even if a**3 is replaced by a*a*a... 484 ENDIF ! aicen > epsi10483 ENDIF 485 484 ! 486 END DO ! ji487 END DO !jj488 END DO !jl485 END DO 486 END DO 487 END DO 489 488 490 489 zzc = rn_pe_rdg * Cp ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation … … 510 509 ! CAN BE REMOVED 511 510 ! 512 IF ( nn_icestr_bvf == 1) THEN511 IF( ln_icestr_bvf ) THEN 513 512 514 513 DO jj = 1, jpj 515 514 DO ji = 1, jpi 516 IF ( bv_i(ji,jj) .GT.0.0 ) THEN515 IF ( bv_i(ji,jj) > 0.0 ) THEN 517 516 zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 ) 518 517 ELSE … … 520 519 ENDIF 521 520 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0))) 522 END DO ! j523 END DO ! i521 END DO 522 END DO 524 523 525 524 ENDIF … … 539 538 DO jj = 2, jpj - 1 540 539 DO ji = 2, jpi - 1 541 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is 542 ! present 540 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present 543 541 zworka(ji,jj) = 4.0 * strength(ji,jj) & 544 542 & + strength(ji-1,jj) * tmask(ji-1,jj,1) & … … 562 560 CALL lbc_lnk( strength, 'T', 1. ) 563 561 564 ENDIF ! ksmooth562 ENDIF 565 563 566 564 !-------------------- … … 579 577 DO jj = 1, jpj - 1 580 578 DO ji = 1, jpi - 1 581 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT.epsi10) THEN ! ice is present579 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present 582 580 numts_rm = 1 ! number of time steps for the running mean 583 IF ( strp1(ji,jj) .GT.0.0 ) numts_rm = numts_rm + 1584 IF ( strp2(ji,jj) .GT.0.0 ) numts_rm = numts_rm + 1581 IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 582 IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 585 583 zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm 586 584 strp2(ji,jj) = strp1(ji,jj) … … 661 659 DO jj = 1, jpj 662 660 DO ji = 1, jpi 663 IF( a_i(ji,jj,jl) .GT.epsi10 ) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)664 ELSE 661 IF( a_i(ji,jj,jl) > epsi10 ) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 662 ELSE ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 665 663 ENDIF 666 664 END DO … … 716 714 ENDIF ! nn_partfun 717 715 718 IF( nn_rafting == 1) THEN ! Ridging and rafting ice participation functions716 IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions 719 717 ! 720 718 DO jl = 1, jpl 721 719 DO jj = 1, jpj 722 720 DO ji = 1, jpi 723 IF ( athorn(ji,jj,jl) .GT.0._wp ) THEN721 IF ( athorn(ji,jj,jl) > 0._wp ) THEN 724 722 !!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time.... 725 723 aridge(ji,jj,jl) = ( TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) … … 727 725 IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp 728 726 aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 ) 729 ENDIF ! athorn730 END DO ! ji731 END DO ! jj732 END DO ! jl733 734 ELSE ! nn_rafting = 0727 ENDIF 728 END DO 729 END DO 730 END DO 731 732 ELSE 735 733 ! 736 734 DO jl = 1, jpl … … 740 738 ENDIF 741 739 742 IF ( nn_rafting == 1) THEN743 744 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT.epsi10 ) THEN740 IF( ln_rafting ) THEN 741 742 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 ) THEN 745 743 DO jl = 1, jpl 746 744 DO jj = 1, jpj 747 745 DO ji = 1, jpi 748 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT.epsi10 ) THEN746 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN 749 747 WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 750 748 WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl … … 792 790 DO ji = 1, jpi 793 791 794 IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT.0.0 ) THEN792 IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN 795 793 hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 796 794 hrmean = MAX(SQRT(rn_hstar*hi), hi*krdgmin) … … 988 986 large_afrft = .false. 989 987 990 !CDIR NODEP991 988 DO ij = 1, icells 992 989 ji = indxi(ij) … … 1115 1112 !-------------------------------------------------------------------- 1116 1113 DO jk = 1, nlay_i 1117 !CDIR NODEP1118 1114 DO ij = 1, icells 1119 1115 ji = indxi(ij) … … 1142 1138 IF( con_i ) THEN 1143 1139 DO jk = 1, nlay_i 1144 !CDIR NODEP1145 1140 DO ij = 1, icells 1146 1141 ji = indxi(ij) … … 1152 1147 1153 1148 IF( large_afrac ) THEN ! there is a bug 1154 !CDIR NODEP1155 1149 DO ij = 1, icells 1156 1150 ji = indxi(ij) … … 1164 1158 ENDIF 1165 1159 IF( large_afrft ) THEN ! there is a bug 1166 !CDIR NODEP1167 1160 DO ij = 1, icells 1168 1161 ji = indxi(ij) … … 1182 1175 DO jl2 = 1, jpl 1183 1176 ! over categories to which ridged ice is transferred 1184 !CDIR NODEP1185 1177 DO ij = 1, icells 1186 1178 ji = indxi(ij) … … 1215 1207 ! Transfer ice energy to category jl2 by ridging 1216 1208 DO jk = 1, nlay_i 1217 !CDIR NODEP1218 1209 DO ij = 1, icells 1219 1210 ji = indxi(ij) … … 1227 1218 DO jl2 = 1, jpl 1228 1219 1229 !CDIR NODEP1230 1220 DO ij = 1, icells 1231 1221 ji = indxi(ij) … … 1242 1232 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + smrft (ji,jj) 1243 1233 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj) 1244 ENDIF ! hraft1234 ENDIF 1245 1235 ! 1246 END DO ! ij1236 END DO 1247 1237 1248 1238 ! Transfer rafted ice energy to category jl2 1249 1239 DO jk = 1, nlay_i 1250 !CDIR NODEP1251 1240 DO ij = 1, icells 1252 1241 ji = indxi(ij) … … 1256 1245 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 1257 1246 ENDIF 1258 END DO ! ij1259 END DO !jk1260 1261 END DO ! jl21247 END DO 1248 END DO 1249 1250 END DO 1262 1251 1263 1252 END DO ! jl1 (deforming categories) … … 1331 1320 !!------------------------------------------------------------------- 1332 1321 INTEGER :: ios ! Local integer output status for namelist read 1333 NAMELIST/namiceitdme/ rn_cs, rn_ pe_rdg, rn_fsnowrdg, rn_fsnowrft, &1334 & rn_gstar, rn_astar, rn_hstar, nn_rafting, rn_hraft, rn_craft, rn_por_rdg, &1322 NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft, & 1323 & rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, & 1335 1324 & nn_partfun 1336 1325 !!------------------------------------------------------------------- … … 1349 1338 WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution ' 1350 1339 WRITE(numout,*)' ~~~~~~~~~~~~~~~' 1351 WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs ', rn_cs 1352 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging rn_pe_rdg ', rn_pe_rdg 1353 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg ', rn_fsnowrdg 1354 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft ', rn_fsnowrft 1355 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar ', rn_gstar 1356 WRITE(numout,*)' Equivalent to G* for an exponential part function rn_astar ', rn_astar 1357 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar ', rn_hstar 1358 WRITE(numout,*)' Rafting of ice sheets or not nn_rafting ', nn_rafting 1359 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft ', rn_hraft 1360 WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft ', rn_craft 1361 WRITE(numout,*)' Initial porosity of ridges rn_por_rdg ', rn_por_rdg 1362 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun ', nn_partfun 1340 WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs 1341 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg 1342 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft 1343 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar 1344 WRITE(numout,*)' Equivalent to G* for an exponential part function rn_astar = ', rn_astar 1345 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar 1346 WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting 1347 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft 1348 WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft 1349 WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg 1350 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun 1363 1351 ENDIF 1364 1352 ! -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r5067 r5070 160 160 zremap_flag(ji,jj) = 0 161 161 ENDIF 162 END DO !ji163 END DO !jj162 END DO 163 END DO 164 164 165 165 !----------------------------------------------------------------------------------------------- … … 201 201 END DO 202 202 203 END DO !jl203 END DO 204 204 205 205 !----------------------------------------------------------------------------------------------- … … 244 244 !----------------------------------------------------------------------------------------------- 245 245 !- 7.1 g(h) for category 1 at start of time step 246 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), & 247 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 246 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 248 247 & hR(:,:,klbnd), zremap_flag ) 249 248 … … 254 253 255 254 !ji 256 IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN255 IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 257 256 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 258 257 ! ji, a_i > epsi10 259 IF (zdh0 .lt. 0.0) THEN !remove area from category 1258 IF( zdh0 < 0.0 ) THEN !remove area from category 1 260 259 ! ji, a_i > epsi10; zdh0 < 0 261 zdh0 = MIN( -zdh0,hi_max(klbnd))260 zdh0 = MIN( -zdh0, hi_max(klbnd) ) 262 261 263 262 !Integrate g(1) from 0 to dh0 to estimate area melted 264 zetamax = MIN( zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd)265 IF (zetamax.gt.0.0) THEN263 zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 264 IF( zetamax > 0.0 ) THEN 266 265 zx1 = zetamax 267 zx2 = 0.5 * zetamax *zetamax266 zx2 = 0.5 * zetamax * zetamax 268 267 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 269 268 ! Constrain new thickness <= ht_i 270 zdamax = a_i(ii,ij,klbnd) * & 271 (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 269 zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! zdamax > 0 272 270 !ice area lost due to melting of thin ice 273 zda0 = MIN( zda0, zdamax)271 zda0 = MIN( zda0, zdamax ) 274 272 275 273 ! Remove area, conserving volume 276 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 277 * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 274 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 278 275 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 279 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) *ht_i(ii,ij,klbnd) ! clem-useless ?280 ENDIF ! zetamax > 0276 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? 277 ENDIF 281 278 ! ji, a_i > epsi10 282 279 283 280 ELSE ! if ice accretion 284 281 ! ji, a_i > epsi10; zdh0 > 0 285 zhbnew(ii,ij,klbnd-1) = MIN( zdh0,hi_max(klbnd))282 zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 286 283 ! zhbnew was 0, and is shifted to the right to account for thin ice 287 284 ! growth in openwater (F0 = f1) … … 295 292 !- 7.3 g(h) for each thickness category 296 293 DO jl = klbnd, kubnd 297 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),&298 g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag)294 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 295 & g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 299 296 END DO 300 297 … … 316 313 ij = nind_j(ji) 317 314 318 IF (zhbnew(ii,ij,jl) .gt.hi_max(jl)) THEN ! transfer from jl to jl+1315 IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 319 316 320 317 ! left and right integration limits in eta space 321 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl)322 zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl)318 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) 319 zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 323 320 zdonor(ii,ij,jl) = jl 324 321 … … 327 324 ! left and right integration limits in eta space 328 325 zvetamin(ji) = 0.0 329 zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1)326 zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1) 330 327 zdonor(ii,ij,jl) = jl + 1 331 328 332 329 ENDIF ! zhbnew(jl) > hi_max(jl) 333 330 334 zetamax = MAX( zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin331 zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 335 332 zetamin = zvetamin(ji) 336 333 337 334 zx1 = zetamax - zetamin 338 zwk1 = zetamin *zetamin339 zwk2 = zetamax *zetamax340 zx2 = 0.5 * ( zwk2 - zwk1)335 zwk1 = zetamin * zetamin 336 zwk2 = zetamax * zetamax 337 zx2 = 0.5 * ( zwk2 - zwk1 ) 341 338 zwk1 = zwk1 * zetamin 342 339 zwk2 = zwk2 * zetamax 343 zx3 = 1.0 /3.0 * (zwk2 - zwk1)340 zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 ) 344 341 nd = zdonor(ii,ij,jl) 345 342 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 346 343 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 347 344 348 END DO ! ji345 END DO 349 346 END DO ! jl klbnd -> kubnd - 1 350 347 … … 365 362 ht_i(ii,ij,1) = rn_himin 366 363 ENDIF 367 END DO !ji364 END DO 368 365 369 366 !!---------------------------------------------------------------------------------------------- … … 401 398 402 399 403 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, & 404 & g0, g1, hL, hR, zremap_flag ) 400 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 405 401 !!------------------------------------------------------------------ 406 402 !! *** ROUTINE lim_itd_fitline *** … … 442 438 ! Change hL or hR if hice falls outside central third of range 443 439 444 zh13 = 1.0 /3.0 * (2.0*hL(ji,jj) + hR(ji,jj))445 zh23 = 1.0 /3.0 * (hL(ji,jj) + 2.0*hR(ji,jj))440 zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 441 zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 446 442 447 443 IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) … … 454 450 zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 455 451 zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 456 g0(ji,jj) = zwk1 * ( 2._wp /3._wp - zwk2 )457 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5)452 g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) 453 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 458 454 ! 459 455 ELSE ! remap_flag = .false. or a_i < epsi10 … … 516 512 517 513 DO jl = klbnd, kubnd 518 zaTsfn(:,:,jl) = a_i(:,:,jl) *t_su(:,:,jl)514 zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 519 515 END DO 520 516 … … 539 535 DO ji = 1, jpi 540 536 541 IF (zdonor(ji,jj,jl) .GT.0) THEN537 IF (zdonor(ji,jj,jl) > 0) THEN 542 538 jl1 = zdonor(ji,jj,jl) 543 539 544 IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 545 IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN 546 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) ) & 547 .OR. & 548 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 549 ) THEN 540 IF (zdaice(ji,jj,jl) < 0.0) THEN 541 IF (zdaice(ji,jj,jl) > -epsi10) THEN 542 IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. & 543 ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 550 544 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 551 545 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 559 553 ENDIF 560 554 561 IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 562 IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN 563 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) ) & 564 .OR. & 565 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 566 ) THEN 555 IF (zdvice(ji,jj,jl) < 0.0) THEN 556 IF (zdvice(ji,jj,jl) > -epsi10 ) THEN 557 IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. & 558 ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 567 559 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 568 560 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 577 569 578 570 ! If daice is close to aicen, set daice = aicen. 579 IF (zdaice(ji,jj,jl) .GT.a_i(ji,jj,jl1) - epsi10 ) THEN580 IF (zdaice(ji,jj,jl) .LT.a_i(ji,jj,jl1)+epsi10) THEN571 IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN 572 IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN 581 573 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 582 574 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 586 578 ENDIF 587 579 588 IF (zdvice(ji,jj,jl) .GT.v_i(ji,jj,jl1)-epsi10) THEN589 IF (zdvice(ji,jj,jl) .LT.v_i(ji,jj,jl1)+epsi10) THEN580 IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN 581 IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN 590 582 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 591 583 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 596 588 597 589 ENDIF ! donor > 0 598 END DO ! i599 END DO ! j600 601 END DO !jl590 END DO 591 END DO 592 593 END DO 602 594 603 595 !------------------------------------------------------------------------------- … … 609 601 DO jj = 1, jpj 610 602 DO ji = 1, jpi 611 IF (zdaice(ji,jj,jl) .GT.0.0 ) THEN ! daice(n) can be < puny603 IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny 612 604 nbrem = nbrem + 1 613 605 nind_i(nbrem) = ji 614 606 nind_j(nbrem) = jj 615 ENDIF ! tmask607 ENDIF 616 608 END DO 617 609 END DO … … 622 614 623 615 jl1 = zdonor(ii,ij,jl) 624 rswitch 625 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * rswitch616 rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 617 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 626 618 IF( jl1 == jl) THEN ; jl2 = jl1+1 627 ELSE 619 ELSE ; jl2 = jl 628 620 ENDIF 629 621 … … 682 674 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 683 675 684 END DO ! ji676 END DO 685 677 686 678 !------------------ … … 689 681 690 682 DO jk = 1, nlay_i 691 !CDIR NODEP692 683 DO ji = 1, nbrem 693 684 ii = nind_i(ji) … … 695 686 696 687 jl1 = zdonor(ii,ij,jl) 697 IF (jl1 .EQ.jl) THEN688 IF (jl1 == jl) THEN 698 689 jl2 = jl+1 699 690 ELSE ! n1 = n+1 … … 704 695 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice 705 696 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 706 END DO ! ji707 END DO ! jk697 END DO 698 END DO 708 699 709 700 END DO ! boundaries, 1 to ncat-1 … … 719 710 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 720 711 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 721 rswitch = 1.0 - MAX( 0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes712 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, -v_s(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 722 713 ELSE 723 714 ht_i(ji,jj,jl) = 0._wp 724 715 t_su(ji,jj,jl) = rtt 725 716 ENDIF 726 END DO ! ji727 END DO ! jj728 END DO ! jl717 END DO 718 END DO 719 END DO 729 720 ! 730 721 CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) … … 836 827 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 ) 837 828 ENDIF 838 END DO ! ji839 END DO ! jj829 END DO 830 END DO 840 831 IF(lk_mpp) CALL mpp_max( zshiftflag ) 841 832 … … 870 861 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 871 862 ENDIF 872 END DO ! ji873 END DO ! jj863 END DO 864 END DO 874 865 875 866 IF(lk_mpp) CALL mpp_max( zshiftflag ) … … 890 881 ! a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 891 882 ! ENDIF 892 ! END DO ! ji893 ! END DO ! jj883 ! END DO 884 ! END DO 894 885 ! clem-change end 895 886 896 END DO ! jl887 END DO 897 888 898 889 !------------------------------------------------------------------------------ -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r5067 r5070 291 291 ! include it later 292 292 293 zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) /e1u(ji,jj)294 zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) /e2v(ji,jj)293 zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 294 zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 295 295 296 296 za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx … … 357 357 & ) * r1_e12t(ji,jj) 358 358 359 zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) /e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &360 & - ( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) /e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &359 zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 360 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 361 361 & ) * r1_e12t(ji,jj) 362 362 363 363 ! 364 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) /e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &365 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) /e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &364 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 365 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 366 366 & ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) ) & 367 367 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) … … 391 391 392 392 !- Calculate Delta on corners 393 zddc = ( ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) /e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &394 & + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) /e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &393 zddc = ( ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 394 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 395 395 & ) * r1_e12f(ji,jj) 396 396 397 zdtc = (- ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) /e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &398 & + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) /e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &397 zdtc = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 398 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 399 399 & ) * r1_e12f(ji,jj) 400 400 … … 420 420 !- contribution of zs1, zs2 and zs12 to zf1 421 421 zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 422 & + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) /e2u(ji,jj) &423 & + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) /e1u(ji,jj) &422 & + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj) & 423 & + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj) & 424 424 & ) * r1_e12u(ji,jj) 425 425 ! contribution of zs1, zs2 and zs12 to zf2 426 426 zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 427 & - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) /e1v(ji,jj) &428 & + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) /e2v(ji,jj) &427 & - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj) & 428 & + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj) & 429 429 & ) * r1_e12v(ji,jj) 430 430 END DO … … 609 609 & ) * r1_e12t(ji,jj) 610 610 611 zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) /e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &612 & -( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) /e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &611 zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 612 & -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 613 613 & ) * r1_e12t(ji,jj) 614 614 ! 615 615 ! SB modif because ocean has no slip boundary condition 616 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) /e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &617 & +( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) /e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &616 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 617 & +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 618 618 & ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) ) & 619 619 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) … … 677 677 DO jj = k_j1+1, k_jpj-1 678 678 DO ji = 2, jpim1 679 IF (zpresh(ji,jj) .GT.1.0) THEN679 IF (zpresh(ji,jj) > 1.0) THEN 680 680 sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 681 681 sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5067 r5070 631 631 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - a_i_1d(ji) + epsi20 ) ) 632 632 ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 633 ! adjust thickness633 ! retrieve total concentration 634 634 at_i_1d(ji) = a_i_1d(ji) 635 635 END IF … … 651 651 !!------------------------------------------------------------------- 652 652 INTEGER :: ios ! Local integer output status for namelist read 653 NAMELIST/namicethd/ rn_hnewice, nn_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, & 654 & rn_himin, parsub, rn_betas, & 655 & rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 653 NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, & 654 & rn_himin, parsub, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 656 655 & nn_monocat 657 656 !!------------------------------------------------------------------- … … 682 681 WRITE(numout,*) 683 682 WRITE(numout,*)' Namelist of ice parameters for ice thermodynamic computation ' 684 WRITE(numout,*)' ice thick. for lateral accretion rn_hnewice 685 WRITE(numout,*)' Frazil ice thickness as a function of wind or not nn_frazil = ', nn_frazil686 WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom rn_maxfrazb 687 WRITE(numout,*)' Thresold relative drift speed for collection of frazil rn_vfrazb 688 WRITE(numout,*)' Squeezing coefficient for collection of frazil rn_Cfrazb 689 WRITE(numout,*)' minimum ice thickness rn_himin 683 WRITE(numout,*)' ice thick. for lateral accretion rn_hnewice = ', rn_hnewice 684 WRITE(numout,*)' Frazil ice thickness as a function of wind or not ln_frazil = ', ln_frazil 685 WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom rn_maxfrazb = ', rn_maxfrazb 686 WRITE(numout,*)' Thresold relative drift speed for collection of frazil rn_vfrazb = ', rn_vfrazb 687 WRITE(numout,*)' Squeezing coefficient for collection of frazil rn_Cfrazb = ', rn_Cfrazb 688 WRITE(numout,*)' minimum ice thickness rn_himin = ', rn_himin 690 689 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 691 690 WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub 692 WRITE(numout,*)' coefficient for ice-lead partition of snowfall rn_betas 693 WRITE(numout,*)' extinction radiation parameter in sea ice rn_kappa_i 691 WRITE(numout,*)' coefficient for ice-lead partition of snowfall rn_betas = ', rn_betas 692 WRITE(numout,*)' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i 694 693 WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nn_conv_dif = ', nn_conv_dif 695 694 WRITE(numout,*)' maximal err. on T for heat diffusion computation rn_terr_dif = ', rn_terr_dif 696 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice nn_ice_thcon 695 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice nn_ice_thcon = ', nn_ice_thcon 697 696 WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i 698 697 WRITE(numout,*)' virtual ITD mono-category parameterizations (1) or not nn_monocat = ', nn_monocat -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5067 r5070 276 276 277 277 ! updates available heat + thickness 278 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )278 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 279 279 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 280 280 … … 573 573 574 574 ENDIF 575 END DO ! ji576 END DO ! jk575 END DO 576 END DO 577 577 578 578 !------------------------------------------- -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5067 r5070 517 517 518 518 DO ji = kideb , kiut 519 IF ( ht_s_1d(ji) .gt.0.0 ) THEN519 IF ( ht_s_1d(ji) > 0.0 ) THEN 520 520 ! 521 521 !------------------------------------------------------------------------------| … … 541 541 ENDIF 542 542 543 IF ( t_su_1d(ji) .LT.rtt ) THEN543 IF ( t_su_1d(ji) < rtt ) THEN 544 544 545 545 !------------------------------------------------------------------------------| … … 587 587 !------------------------------------------------------------------------------| 588 588 ! 589 IF (t_su_1d(ji) .LT.rtt) THEN589 IF (t_su_1d(ji) < rtt) THEN 590 590 ! 591 591 !------------------------------------------------------------------------------| … … 705 705 DO ji = kideb , kiut 706 706 ! snow temperatures 707 IF (ht_s_1d(ji) .GT.0._wp) &707 IF (ht_s_1d(ji) > 0._wp) & 708 708 t_s_1d(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 709 709 * t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r5067 r5070 156 156 hicol(:,:) = rn_hnewice 157 157 158 IF( nn_frazil == 1) THEN158 IF( ln_frazil ) THEN 159 159 160 160 !-------------------- … … 221 221 iterate_frazil = .true. 222 222 223 DO WHILE ( iter .LT.100 .AND. iterate_frazil )223 DO WHILE ( iter < 100 .AND. iterate_frazil ) 224 224 zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) & 225 225 - hicol(ji,jj) * zhicrit * ztwogp * zvrel2 … … 321 321 zh_newice(ji) = rn_hnewice 322 322 END DO 323 IF( nn_frazil == 1) zh_newice(1:nbpac) = hicol_1d(1:nbpac)323 IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 324 324 325 325 !---------------------- -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
r5067 r5070 150 150 WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity ' 151 151 WRITE(numout,*) '~~~~~~~~~~~~~~~~' 152 WRITE(numout,*) ' switch for salinity nn_icesal :', nn_icesal153 WRITE(numout,*) ' bulk salinity value if nn_icesal = 1 :', rn_icesal154 WRITE(numout,*) ' restoring salinity for GD :', rn_sal_gd155 WRITE(numout,*) ' restoring time for GD :', rn_time_gd156 WRITE(numout,*) ' restoring salinity for flushing :', rn_sal_fl157 WRITE(numout,*) ' restoring time for flushing :', rn_time_fl158 WRITE(numout,*) ' Maximum tolerated ice salinity :', rn_simax159 WRITE(numout,*) ' Minimum tolerated ice salinity :', rn_simin152 WRITE(numout,*) ' switch for salinity nn_icesal = ', nn_icesal 153 WRITE(numout,*) ' bulk salinity value if nn_icesal = 1 = ', rn_icesal 154 WRITE(numout,*) ' restoring salinity for GD = ', rn_sal_gd 155 WRITE(numout,*) ' restoring time for GD = ', rn_time_gd 156 WRITE(numout,*) ' restoring salinity for flushing = ', rn_sal_fl 157 WRITE(numout,*) ' restoring time for flushing = ', rn_time_fl 158 WRITE(numout,*) ' Maximum tolerated ice salinity = ', rn_simax 159 WRITE(numout,*) ' Minimum tolerated ice salinity = ', rn_simin 160 160 ENDIF 161 161 ! -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r5067 r5070 106 106 ! mass and salt flux init 107 107 zviold(:,:,:) = v_i(:,:,:) 108 zvsold(:,:,:) = v_s(:,:,:) 108 109 zeiold(:,:) = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 109 110 zesold(:,:) = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) … … 131 132 132 133 ! If ice drift field is too fast, use an appropriate time step for advection. 133 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice /e1u(:,:) ) ! CFL test for stability134 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice /e2v(:,:) ) )134 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) ! CFL test for stability 135 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 135 136 IF(lk_mpp ) CALL mpp_max( zcfl ) 136 137 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r5067 r5070 26 26 REAL(wp), PUBLIC :: rn_hnewice !: thickness for new ice formation (m) 27 27 28 INTEGER , PUBLIC :: nn_frazil !: use of frazil ice collection as function of wind (1) or not (0)28 LOGICAL , PUBLIC :: ln_frazil !: use of frazil ice collection as function of wind (T) or not (F) 29 29 30 30 !!----------------------------- -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r5067 r5070 369 369 WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 370 370 WRITE(numout,*) ' ~~~~~~' 371 WRITE(numout,*) ' number of ice categories 372 WRITE(numout,*) ' number of ice layers 373 WRITE(numout,*) ' number of snow layers 371 WRITE(numout,*) ' number of ice categories = ', jpl 372 WRITE(numout,*) ' number of ice layers = ', nlay_i 373 WRITE(numout,*) ' number of snow layers = ', nlay_s 374 374 WRITE(numout,*) ' switch for ice dynamics (1) or not (0) ln_limdyn = ', ln_limdyn 375 375 WRITE(numout,*) ' maximum ice concentration = ', rn_amax
Note: See TracChangeset
for help on using the changeset viewer.