Changeset 8313 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
- Timestamp:
- 2017-07-10T20:24:21+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 11 edited
-
ice.F90 (modified) (5 diffs)
-
limctl.F90 (modified) (1 diff)
-
limdyn.F90 (modified) (3 diffs)
-
limrhg.F90 (modified) (4 diffs)
-
limsbc.F90 (modified) (4 diffs)
-
limthd.F90 (modified) (7 diffs)
-
limthd_dh.F90 (modified) (1 diff)
-
limthd_dif.F90 (modified) (14 diffs)
-
limthd_lac.F90 (modified) (1 diff)
-
limwri.F90 (modified) (1 diff)
-
thd_ice.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r8312 r8313 196 196 197 197 ! !!** ice-thickness distribution namelist (namiceitd) ** 198 INTEGER , PUBLIC :: nn_catbnd !: categories distribution following: tanh function (1), or h^(-alpha) function (2)199 198 REAL(wp), PUBLIC :: rn_himean !: mean thickness of the domain (used to compute the distribution, nn_itdshp = 2 only) 200 199 … … 210 209 LOGICAL , PUBLIC :: ln_icestr_bvf !: use brine volume to diminish ice strength 211 210 ! -- limdyn & limrhg -- ! 211 REAL(wp), PUBLIC :: rn_ishlat !: lateral boundary condition for sea-ice 212 212 REAL(wp), PUBLIC :: rn_cio !: drag coefficient for oceanic stress 213 213 REAL(wp), PUBLIC :: rn_creepl !: creep limit : has to be under 1.0e-9 … … 223 223 ! -- limthd_dif -- ! 224 224 REAL(wp), PUBLIC :: rn_kappa_i !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 225 REAL(wp), PUBLIC :: nn_conv_dif !: maximal number of iterations for heat diffusion226 REAL(wp), PUBLIC :: rn_terr_dif !: maximal tolerated error (C) for heat diffusion227 225 INTEGER , PUBLIC :: nn_ice_thcon !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007) 228 LOGICAL , PUBLIC :: ln_ it_qnsice !: iterate surface flux with changing surface temperatureor not (F)229 INTEGER , PUBLIC :: nn_monocat !: virtual ITD mono-category parameterizations (1 ) or not (0)226 LOGICAL , PUBLIC :: ln_dqnsice !: change non-solar surface flux with changing surface temperature (T) or not (F) 227 INTEGER , PUBLIC :: nn_monocat !: virtual ITD mono-category parameterizations (1-4) or not (0) 230 228 REAL(wp), PUBLIC :: rn_cdsn !: thermal conductivity of the snow [W/m/K] 231 229 ! -- limthd_dh -- ! … … 310 308 ! 311 309 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_bo !: Sea-Ice bottom temperature [Kelvin] 312 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: frld !: Leads fraction = 1 - ice fraction313 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: pfrld !: Leads fraction at previous time314 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: phicif !: Old ice thickness315 310 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qlead !: heat balance of the lead (or of the open ocean) 316 311 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fhtur !: net downward heat flux from the ice to the ocean … … 534 529 535 530 ii = ii + 1 536 ALLOCATE( t_bo (jpi,jpj) , frld (jpi,jpj) , pfrld (jpi,jpj) , phicif (jpi,jpj) ,&531 ALLOCATE( t_bo (jpi,jpj) , & 537 532 & wfx_snw(jpi,jpj) , wfx_snw_dyn(jpi,jpj) , wfx_snw_sum(jpi,jpj) , wfx_snw_sub(jpi,jpj) , & 538 533 & wfx_ice(jpi,jpj) , wfx_sub (jpi,jpj) , wfx_ice_sub(jpi,jpj) , wfx_lam (jpi,jpj) , & -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
r7646 r8313 385 385 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ ' 386 386 WRITE(numout,*) ' - Heat fluxes in and out the ice ***' 387 WRITE(numout,*) ' qsr_ini : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) )388 WRITE(numout,*) ' qns_ini : ', pfrld(ji,jj) * qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) )387 WRITE(numout,*) ' qsr_ini : ', (1._wp-at_i_b(ji,jj)) * qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) ) 388 WRITE(numout,*) ' qns_ini : ', (1._wp-at_i_b(ji,jj)) * qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) ) 389 389 WRITE(numout,*) 390 390 WRITE(numout,*) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r7753 r8313 109 109 !!------------------------------------------------------------------- 110 110 INTEGER :: ios ! Local integer output status for namelist read 111 NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord, & 112 & nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, & 111 NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord, & 112 & nn_icestr, rn_pe_rdg, rn_pstar, rn_crhg, ln_icestr_bvf, & 113 & rn_ishlat, rn_cio, rn_creepl, rn_ecc, & 113 114 & nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr, rn_lfrelax 114 115 !!------------------------------------------------------------------- … … 130 131 WRITE(numout,*)' choose the advection scheme (-1=Prather, 0=Ulimate-Macho) nn_limadv = ', nn_limadv 131 132 WRITE(numout,*)' choose the order of the scheme (if ultimate) nn_limadv_ord = ', nn_limadv_ord 132 ! lim rhg133 ! limitd_me 133 134 WRITE(numout,*)' ice strength parameterization (0=Hibler 1=Rothrock) nn_icestr = ', nn_icestr 134 WRITE(numout,*)' Including brine volume in ice strength comp. ln_icestr_bvf = ', ln_icestr_bvf135 135 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging rn_pe_rdg = ', rn_pe_rdg 136 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio137 136 WRITE(numout,*) ' first bulk-rheology parameter rn_pstar = ', rn_pstar 138 137 WRITE(numout,*) ' second bulk-rhelogy parameter rn_crhg = ', rn_crhg 138 WRITE(numout,*)' Including brine volume in ice strength comp. ln_icestr_bvf = ', ln_icestr_bvf 139 ! limrhg 140 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat 141 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 139 142 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 140 143 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc … … 142 145 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 143 146 WRITE(numout,*) ' Landfast: param (T or F) ln_landfast = ', ln_landfast 144 WRITE(numout,*) ' Landfast: fraction of ocean depth that ice must reach rn_gamma = ', rn_gamma 145 WRITE(numout,*) ' Landfast: maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr 146 WRITE(numout,*) ' Landfast: relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax 147 WRITE(numout,*) ' T: fraction of ocean depth that ice must reach rn_gamma = ', rn_gamma 148 WRITE(numout,*) ' T: maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr 149 WRITE(numout,*) ' T: relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax 150 ENDIF 151 ! 152 IF ( rn_ishlat == 0. ) THEN ; IF(lwp) WRITE(numout,*) ' ice lateral free-slip ' 153 ELSEIF ( rn_ishlat == 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ice lateral no-slip ' 154 ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ice lateral partial-slip ' 155 ELSEIF ( 2. < rn_ishlat ) THEN ; IF(lwp) WRITE(numout,*) ' ice lateral strong-slip ' 147 156 ENDIF 148 157 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r8306 r8313 150 150 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 151 151 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity 152 REAL(wp), PARAMETER :: zshlat = 2._wp ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)153 152 !!------------------------------------------------------------------- 154 153 … … 182 181 DO ji = fs_2, fs_jpim1 ! vector opt. 183 182 IF( zfmask(ji,jj) == 0._wp ) THEN 184 zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )183 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 185 184 ENDIF 186 185 END DO … … 188 187 DO jj = 2, jpjm1 189 188 IF( zfmask(1,jj) == 0._wp ) THEN 190 zfmask(1 ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )189 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 191 190 ENDIF 192 191 IF( zfmask(jpi,jj) == 0._wp ) THEN 193 zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )192 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 194 193 ENDIF 195 194 END DO 196 195 DO ji = 2, jpim1 197 196 IF( zfmask(ji,1) == 0._wp ) THEN 198 zfmask(ji,1 ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )197 zfmask(ji,1 ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 199 198 ENDIF 200 199 IF( zfmask(ji,jpj) == 0._wp ) THEN 201 zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )200 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 202 201 ENDIF 203 202 END DO … … 251 250 !------------------------------------------------------------------------------! 252 251 253 IF( nn_ice_embd == 2) THEN !== embedded sea ice: compute representative ice top surface ==!252 IF( ln_ice_embd ) THEN !== embedded sea ice: compute representative ice top surface ==! 254 253 ! 255 254 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r8266 r8313 112 112 ! --- case we bypass ice thermodynamics --- ! 113 113 IF( .NOT. ln_limthd ) THEN ! we suppose ice is impermeable => ocean is isolated from atmosphere 114 hfx_in (:,:) = pfrld(:,:) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)115 hfx_out (:,:) = pfrld(:,:) * qns_oce(:,:) + qemp_oce(:,:)114 hfx_in (:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:) 115 hfx_out (:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) 116 116 ftr_ice (:,:,:) = 0._wp 117 117 emp_ice (:,:) = 0._wp … … 195 195 & + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) + sfx_sub(:,:) + sfx_lam(:,:) 196 196 197 !-------------------------------------------------------------! 198 ! mass of snow and ice per unit area for embedded sea-ice ! 199 !-------------------------------------------------------------! 200 IF( nn_ice_embd /= 0 ) THEN 201 ! save mass from the previous ice time step 202 snwice_mass_b(:,:) = snwice_mass(:,:) 203 ! new mass per unit area 204 snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 205 ! time evolution of snow+ice mass 206 snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice 207 ENDIF 197 !----------------------------------------! 198 ! mass of snow and ice per unit area ! 199 !----------------------------------------! 200 ! save mass from the previous ice time step 201 snwice_mass_b(:,:) = snwice_mass(:,:) 202 ! new mass per unit area 203 snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 204 ! time evolution of snow+ice mass 205 snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice 208 206 209 207 !-----------------------------------------------! … … 243 241 !! tmod_io = rhoco * | U_ice-U_oce | 244 242 !! - update the modulus of stress at ocean surface 245 !! taum = frld * taum + (1-frld)* tmod_io * | U_ice-U_oce |243 !! taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce | 246 244 !! * at each ocean time step (every kt): 247 245 !! compute linearized ice-ocean stresses as … … 336 334 ! 337 335 IF( .NOT. ln_rstart ) THEN 338 ! ! embedded sea ice 339 IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass 340 snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 341 snwice_mass_b(:,:) = snwice_mass(:,:) 342 ELSE 343 snwice_mass (:,:) = 0._wp ! no mass exchanges 344 snwice_mass_b(:,:) = 0._wp ! no mass exchanges 345 ENDIF 346 IF( nn_ice_embd == 2 ) THEN ! full embedment (case 2) deplete the initial ssh below sea-ice area 336 ! 337 snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) ! snow+ice mass 338 snwice_mass_b(:,:) = snwice_mass(:,:) 339 ! 340 IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area 347 341 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 348 342 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r8239 r8313 172 172 ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 173 173 zqld = tmask(ji,jj,1) * rdt_ice * & 174 & ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 174 & ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) + & 175 & ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 175 176 176 177 ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! … … 210 211 IF( .NOT. ln_limdO ) qlead(:,:) = 0._wp 211 212 ! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere 212 IF( .NOT. ln_limdH ) hfx_in(:,:) = pfrld(:,:) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)213 IF( .NOT. ln_limdH ) hfx_in(:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:) 213 214 IF( .NOT. ln_limdH ) fhtur (:,:) = 0._wp ; fhld (:,:) = 0._wp 214 215 … … 221 222 DO jj = 1, jpj 222 223 DO ji = 1, jpi 223 hfx_out(ji,jj) = pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) & ! Non solar heat flux received by the ocean224 & - qlead(ji,jj) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation225 & - at_i(ji,jj) * fhtur(ji,jj) & ! heat flux taken by turbulence226 & - at_i(ji,jj) * fhld(ji,jj) ! heat flux taken during bottom growth/melt227 ! (fhld should be 0 while bott growth)224 hfx_out(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) & ! Non solar heat flux received by the ocean 225 & - qlead(ji,jj) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation 226 & - at_i(ji,jj) * fhtur(ji,jj) & ! heat flux taken by turbulence 227 & - at_i(ji,jj) * fhld(ji,jj) ! heat flux taken during bottom growth/melt 228 ! (fhld should be 0 while bott growth) 228 229 END DO 229 230 END DO … … 590 591 !!------------------------------------------------------------------- 591 592 INTEGER :: ios ! Local integer output status for namelist read 592 NAMELIST/namicethd/ rn_kappa_i, nn_ conv_dif, rn_terr_dif, nn_ice_thcon, ln_it_qnsice, nn_monocat, rn_cdsn,&593 & ln_limdH, rn_betas, 594 & ln_limdA, rn_beta, rn_dmin, 593 NAMELIST/namicethd/ rn_kappa_i, nn_ice_thcon, ln_dqnsice, rn_cdsn, & 594 & ln_limdH, rn_betas, & 595 & ln_limdA, rn_beta, rn_dmin, & 595 596 & ln_limdO, rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, rn_himin 596 597 !!------------------------------------------------------------------- … … 605 606 IF(lwm) WRITE ( numoni, namicethd ) 606 607 ! 607 IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN608 nn_monocat = 0609 IF(lwp) WRITE(numout,*)610 IF(lwp) WRITE(numout,*) ' nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'611 ENDIF612 608 ! 613 609 IF(lwp) THEN ! control print … … 616 612 WRITE(numout,*)' -- limthd_dif --' 617 613 WRITE(numout,*)' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i 618 WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nn_conv_dif = ', nn_conv_dif619 WRITE(numout,*)' maximal err. on T for heat diffusion computation rn_terr_dif = ', rn_terr_dif620 614 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice nn_ice_thcon = ', nn_ice_thcon 621 WRITE(numout,*)' iterate the surface non-solar flux (T) or not (F) ln_it_qnsice = ', ln_it_qnsice 622 WRITE(numout,*)' virtual ITD mono-category parameterizations (1) or not nn_monocat = ', nn_monocat 615 WRITE(numout,*)' change the surface non-solar flux with Tsu or not ln_dqnsice = ', ln_dqnsice 623 616 WRITE(numout,*)' thermal conductivity of the snow rn_cdsn = ', rn_cdsn 624 617 WRITE(numout,*)' -- limthd_dh --' … … 641 634 ENDIF 642 635 ! 636 IF ( rn_hnewice < rn_himin ) CALL ctl_stop( 'STOP', 'lim_thd_init : rn_hnewice should be >= rn_himin' ) 637 ! 643 638 END SUBROUTINE lim_thd_init 644 639 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r8239 r8313 702 702 !!-------------------------------------------------------------------------- 703 703 SUBROUTINE lim_thd_snwblow_2d( pin, pout ) 704 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pin ! previous fraction lead ( pfrld or (1. - a_i_b))704 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pin ! previous fraction lead ( 1. - a_i_b ) 705 705 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout 706 706 pout = ( 1._wp - ( pin )**rn_betas ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r8239 r8313 92 92 93 93 !! * Local variables 94 INTEGER :: ji ! spatial loop index 95 INTEGER :: ii, ij ! temporary dummy loop index 96 INTEGER :: numeq ! current reference number of equation 97 INTEGER :: jk ! vertical dummy loop index 98 INTEGER :: nconv ! number of iterations in iterative procedure 94 INTEGER :: ji ! spatial loop index 95 INTEGER :: ii, ij ! temporary dummy loop index 96 INTEGER :: numeq ! current reference number of equation 97 INTEGER :: jk ! vertical dummy loop index 99 98 INTEGER :: minnumeqmin, maxnumeqmax 99 INTEGER :: iconv ! number of iterations in iterative procedure 100 INTEGER :: iconv_max = 50 ! max number of iterations in iterative procedure 100 101 101 102 INTEGER, POINTER, DIMENSION(:) :: numeqmin ! reference number of top equation … … 109 110 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 110 111 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered as 0°C 111 REAL(wp) :: ztmelt_i ! ice melting temperature 112 REAL(wp) :: zerritmax ! current maximal error on temperature 112 REAL(wp) :: ztmelt_i ! ice melting temperature 113 113 REAL(wp) :: zhsu 114 REAL(wp) :: zdti_max ! current maximal error on temperature 115 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 114 116 115 117 REAL(wp), POINTER, DIMENSION(:) :: isnow ! switch for presence (1) or absence (0) of snow … … 122 124 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 123 125 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function 124 REAL(wp), POINTER, DIMENSION(:) :: z errit! current error on temperature126 REAL(wp), POINTER, DIMENSION(:) :: zdti ! current error on temperature 125 127 REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4) 126 128 REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice … … 168 170 CALL wrk_alloc( jpij, numeqmin, numeqmax ) 169 171 CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 170 CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, z errit, zdifcase, zftrice, zihic, zghe )172 CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zdti, zdifcase, zftrice, zihic, zghe ) 171 173 CALL wrk_alloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0 ) 172 174 CALL wrk_alloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0 ) … … 290 292 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter 291 293 t_su_1d(ji) = MIN( t_su_1d(ji), rt0 - ztsu_err ) ! necessary 292 z errit(ji) = 1000._wp ! initial value of error294 zdti (ji) = 1000._wp ! initial value of error 293 295 END DO 294 296 … … 305 307 END DO 306 308 307 nconv= 0 ! number of iterations308 z erritmax = 1000._wp ! maximal value of error on all points309 310 DO WHILE ( z erritmax > rn_terr_dif .AND. nconv < nn_conv_dif)311 ! 312 nconv = nconv + 1309 iconv = 0 ! number of iterations 310 zdti_max = 1000._wp ! maximal value of error on all points 311 312 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) 313 ! 314 iconv = iconv + 1 313 315 ! 314 316 !------------------------------------------------------------------------------| … … 452 454 !------------------------------------------------------------------------------| 453 455 ! 454 IF ( ln_ it_qnsice ) THEN456 IF ( ln_dqnsice ) THEN 455 457 DO ji = kideb , kiut 456 458 ! update of the non solar flux according to the update in T_su … … 703 705 ! 704 706 ! check that nowhere it has started to melt 705 ! z errit(ji) is a measure of error, it has to be under terr_dif707 ! zdti(ji) is a measure of error, it has to be under zdti_bnd 706 708 DO ji = kideb , kiut 707 709 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , 190._wp ) 708 z errit(ji)= ABS( t_su_1d(ji) - ztsubit(ji) )710 zdti (ji) = ABS( t_su_1d(ji) - ztsubit(ji) ) 709 711 END DO 710 712 … … 712 714 DO ji = kideb , kiut 713 715 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), 190._wp ) 714 z errit(ji) = MAX( zerrit(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )716 zdti (ji) = MAX( zdti(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) ) 715 717 END DO 716 718 END DO … … 720 722 ztmelt_i = -tmut * s_i_1d(ji,jk) + rt0 721 723 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp ) 722 z errit(ji) = MAX( zerrit(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) )724 zdti (ji) = MAX( zdti(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) ) 723 725 END DO 724 726 END DO … … 726 728 ! Compute spatial maximum over all errors 727 729 ! note that this could be optimized substantially by iterating only the non-converging points 728 z erritmax = 0._wp730 zdti_max = 0._wp 729 731 DO ji = kideb, kiut 730 z erritmax = MAX( zerritmax, zerrit(ji) )731 END DO 732 IF( lk_mpp ) CALL mpp_max( z erritmax, kcom=ncomm_ice )732 zdti_max = MAX( zdti_max, zdti(ji) ) 733 END DO 734 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 733 735 734 736 END DO ! End of the do while iterative procedure … … 745 747 746 748 IF( ln_limctl .AND. lwp ) THEN 747 WRITE(numout,*) ' z erritmax : ', zerritmax748 WRITE(numout,*) ' nconv : ', nconv749 WRITE(numout,*) ' zdti_max : ', zdti_max 750 WRITE(numout,*) ' iconv : ', iconv 749 751 ENDIF 750 752 … … 772 774 773 775 ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 774 IF ( ln_ it_qnsice ) THEN776 IF ( ln_dqnsice ) THEN 775 777 DO ji = kideb, kiut 776 778 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) … … 810 812 CALL wrk_dealloc( jpij, numeqmin, numeqmax ) 811 813 CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 812 CALL wrk_dealloc( jpij, zf, dzf, zqns_ice_b, z errit, zdifcase, zftrice, zihic, zghe )814 CALL wrk_dealloc( jpij, zf, dzf, zqns_ice_b, zdti, zdifcase, zftrice, zihic, zghe ) 813 815 CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 814 816 CALL wrk_dealloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r7753 r8313 67 67 !! - Updating ice internal temperature 68 68 !! - Computation of variation of ice volume and mass 69 !! - Computation of frldbafter lateral accretion and69 !! - Computation of a_i after lateral accretion and 70 70 !! update ht_s_1d, ht_i_1d and tbif_1d(:,:) 71 71 !!------------------------------------------------------------------------ -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r8294 r8313 109 109 !---------------------------------------- 110 110 ! fluxes 111 ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 112 IF( iom_use('qsr_oce') ) CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) ) ! solar flux at ocean surface 113 IF( iom_use('qns_oce') ) CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) ) ! non-solar flux at ocean surface 111 IF( iom_use('qsr_oce') ) CALL iom_put( "qsr_oce" , qsr_oce(:,:) * ( 1._wp - at_i_b(:,:) ) ) ! solar flux at ocean surface 112 IF( iom_use('qns_oce') ) CALL iom_put( "qns_oce" , qns_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) ) ! non-solar flux at ocean surface 114 113 IF( iom_use('qsr_ice') ) CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux at ice surface 115 114 IF( iom_use('qns_ice') ) CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 116 115 IF( iom_use('qtr_ice') ) CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux transmitted thru ice 117 IF( iom_use('qt_oce' ) ) CALL iom_put( "qt_oce" , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )116 IF( iom_use('qt_oce' ) ) CALL iom_put( "qt_oce" , ( qsr_oce(:,:) + qns_oce(:,:) ) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) ) 118 117 IF( iom_use('qt_ice' ) ) CALL iom_put( "qt_ice" , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) & 119 118 & * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r8240 r8313 82 82 83 83 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sprecip_1d !: <==> the 2D sprecip 84 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: frld_1d !: <==> the 2D frld85 84 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: at_i_1d !: <==> the 2D at_i 86 85 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhtur_1d !: <==> the 2D fhtur … … 153 152 ! 154 153 ii = ii + 1 155 ALLOCATE( sprecip_1d (jpij) , frld_1d (jpij) ,at_i_1d (jpij) , &154 ALLOCATE( sprecip_1d (jpij) , at_i_1d (jpij) , & 156 155 & fhtur_1d (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) , wfx_snw_sum_1d(jpij) , & 157 156 & fhld_1d (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) , &
Note: See TracChangeset
for help on using the changeset viewer.
