Changeset 2715 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
- Timestamp:
- 2011-03-30T17:58:35+02:00 (13 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r2591 r2715 5 5 !! computation of surface and inner T 6 6 !!====================================================================== 7 !! History : LIM ! 02-2003 (M. Vancoppenolle) original 1D code 8 !! ! 06-2005 (M. Vancoppenolle) 3d version 9 !! ! 11-2006 (X Fettweis) Vectorization by Xavier 10 !! ! 04-2007 (M. Vancoppenolle) Energy conservation 11 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 7 12 !!---------------------------------------------------------------------- 8 13 #if defined key_lim3 … … 12 17 USE par_oce ! ocean parameters 13 18 USE phycst ! physical constants (ocean directory) 14 USE thd_ice15 USE in_out_manager16 USE ice17 USE par_ice18 USE lib_mpp 19 USE ice ! LIM-3 variables 20 USE par_ice ! LIM-3 parameters 21 USE thd_ice ! LIM-3: thermodynamics 22 USE in_out_manager ! I/O manager 23 USE lib_mpp ! MPP library 19 24 20 25 IMPLICIT NONE … … 23 28 PUBLIC lim_thd_dif ! called by lim_thd 24 29 25 REAL(wp) :: & ! constant values 26 epsi20 = 1e-20 , & 27 epsi13 = 1e-13 , & 28 zzero = 0.e0 , & 29 zone = 1.e0 30 REAL(wp) :: epsi20 = 1e-20 ! constant values 31 REAL(wp) :: epsi13 = 1e-13 ! constant values 30 32 31 33 !!---------------------------------------------------------------------- 32 !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)34 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 33 35 !! $Id$ 34 36 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 95 97 !! * Local variables 96 98 INTEGER :: ji, & ! spatial loop index 97 zji, zjj, & ! temporary dummy loop index99 ii, ij, & ! temporary dummy loop index 98 100 numeq, & ! current reference number of equation 99 101 layer, & ! vertical dummy loop index 100 102 nconv, & ! number of iterations in iterative procedure 101 minnumeqmin, & ! 102 maxnumeqmax 103 minnumeqmin, maxnumeqmax 103 104 104 105 INTEGER , DIMENSION(kiut) :: & … … 137 138 zdiagbis 138 139 139 REAL(wp) , DIMENSION(kiut,jkmax+2,3) :: & 140 ztrid ! tridiagonal system terms 140 REAL(wp) , DIMENSION(kiut,jkmax+2,3) :: ztrid ! tridiagonal system terms 141 141 142 142 REAL(wp), DIMENSION(kiut) :: & 143 143 ztfs , & ! ice melting point 144 ztsuold , & ! old surface temperature (before the iterative 145 ! procedure ) 144 ztsuold , & ! old surface temperature (before the iterative procedure ) 146 145 ztsuoldit, & ! surface temperature at previous iteration 147 146 zh_i , & !ice layer thickness … … 152 151 153 152 REAL(wp) :: & ! constant values 154 zeps = 1.0e-10, & ! 155 zg1s = 2.0, & !: for the tridiagonal system 156 zg1 = 2.0, & 157 zgamma = 18009.0, & !: for specific heat 158 zbeta = 0.117, & !: for thermal conductivity (could be 0.13) 159 zraext_s = 1.0e08, & !: extinction coefficient of radiation in the snow 160 zkimin = 0.10 , & !: minimum ice thermal conductivity 161 zht_smin = 1.0e-4 !: minimum snow depth 162 163 REAL(wp) :: & ! local variables 164 ztmelt_i, & ! ice melting temperature 165 zerritmax ! current maximal error on temperature 166 167 REAL(wp), DIMENSION(kiut) :: & 168 zerrit, & ! current error on temperature 169 zdifcase, & ! case of the equation resolution (1->4) 170 zftrice, & ! solar radiation transmitted through the ice 171 zihic, zhsu 153 zeps = 1.e-10_wp, & ! 154 zg1s = 2._wp, & !: for the tridiagonal system 155 zg1 = 2._wp, & 156 zgamma = 18009._wp, & !: for specific heat 157 zbeta = 0.117_wp, & !: for thermal conductivity (could be 0.13) 158 zraext_s = 1.e+8_wp, & !: extinction coefficient of radiation in the snow 159 zkimin = 0.10_wp , & !: minimum ice thermal conductivity 160 zht_smin = 1.e-4_wp !: minimum snow depth 161 162 REAL(wp) :: ztmelt_i ! ice melting temperature 163 REAL(wp) :: zerritmax ! current maximal error on temperature 164 REAL(wp), DIMENSION(kiut) :: zerrit ! current error on temperature 165 REAL(wp), DIMENSION(kiut) :: zdifcase ! case of the equation resolution (1->4) 166 REAL(wp), DIMENSION(kiut) :: zftrice ! solar radiation transmitted through the ice 167 REAL(wp), DIMENSION(kiut) :: zihic, zhsu 172 168 !!------------------------------------------------------------------ 173 169 ! … … 178 174 DO ji = kideb , kiut 179 175 ! is there snow or not 180 isnow(ji)= INT ( 1.0 - MAX( 0.0 , SIGN (1.0, - ht_s_b(ji) ) ))176 isnow(ji)= INT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) ) 181 177 ! surface temperature of fusion 178 !!gm ??? ztfs(ji) = rtt !!!???? 182 179 ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt 183 180 ! layer thickness 184 zh_i(ji) 185 zh_s(ji) 181 zh_i(ji) = ht_i_b(ji) / nlay_i 182 zh_s(ji) = ht_s_b(ji) / nlay_s 186 183 END DO 187 184 … … 190 187 !-------------------- 191 188 192 z_s(:,0) = 0.0 ! vert. coord. of the up. lim. of the 1st snow layer 193 z_i(:,0) = 0.0 ! vert. coord. of the up. lim. of the 1st ice layer 194 195 DO layer = 1, nlay_s 196 DO ji = kideb , kiut 197 ! vert. coord of the up. lim. of the layer-th snow layer 198 z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s 199 END DO 200 END DO 201 202 DO layer = 1, nlay_i 203 DO ji = kideb , kiut 204 ! vert. coord of the up. lim. of the layer-th ice layer 205 z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i 189 z_s(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st snow layer 190 z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer 191 192 DO layer = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 193 DO ji = kideb , kiut 194 z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s 195 END DO 196 END DO 197 198 DO layer = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 199 DO ji = kideb , kiut 200 z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i 206 201 END DO 207 202 END DO … … 224 219 DO ji = kideb , kiut 225 220 ! switches 226 isnow(ji) = INT ( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ))221 isnow(ji) = INT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ) ) 227 222 ! hs > 0, isnow = 1 228 zhsu(ji) = hnzst !threshold for the computation of i0 229 zihic(ji) = MAX( zzero , 1.0 - ( ht_i_b(ji) / zhsu(ji) ) ) 230 231 i0(ji) = ( 1.0 - isnow(ji) ) * & 232 ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 223 zhsu (ji) = hnzst ! threshold for the computation of i0 224 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) ) 225 226 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 233 227 !fr1_i0_1d = i0 for a thin ice surface 234 228 !fr1_i0_2d = i0 for a thick ice surface … … 244 238 !------------------------------------------------------- 245 239 DO ji = kideb , kiut 246 247 ! Shortwave radiation absorbed at surface 248 zfsw(ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) 249 250 ! Solar radiation transmitted below the surface layer 251 zftrice(ji)= qsr_ice_1d(ji) * i0(ji) 252 253 ! derivative of incoming nonsolar flux 254 dzf(ji) = dqns_ice_1d(ji) 255 240 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 241 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 242 dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 256 243 END DO 257 244 … … 260 247 !--------------------------------------------------------- 261 248 262 DO ji = kideb , kiut 263 ! Initialization 264 zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow 265 END DO 266 267 ! Radiation through snow 268 DO layer = 1, nlay_s 269 DO ji = kideb , kiut 270 ! radiation transmitted below the layer-th snow layer 271 zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP ( - zraext_s * ( MAX ( 0.0 , & 272 z_s(ji,layer) ) ) ) 273 ! radiation absorbed by the layer-th snow layer 249 DO ji = kideb, kiut ! snow initialization 250 zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow 251 END DO 252 253 DO layer = 1, nlay_s ! Radiation through snow 254 DO ji = kideb, kiut 255 ! ! radiation transmitted below the layer-th snow layer 256 zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) ) 257 ! ! radiation absorbed by the layer-th snow layer 274 258 zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer) 275 259 END DO 276 260 END DO 277 261 278 ! Radiation through ice 279 DO ji = kideb , kiut 280 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + & 281 zftrice(ji) * ( 1 - isnow(ji) ) 282 END DO 283 284 DO layer = 1, nlay_i 285 DO ji = kideb , kiut 286 ! radiation transmitted below the layer-th ice layer 287 zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP ( - kappa_i * ( MAX ( 0.0 , & 288 z_i(ji,layer) ) ) ) 289 ! radiation absorbed by the layer-th ice layer 262 DO ji = kideb, kiut ! ice initialization 263 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 264 END DO 265 266 DO layer = 1, nlay_i ! Radiation through ice 267 DO ji = kideb, kiut 268 ! ! radiation transmitted below the layer-th ice layer 269 zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) ) 270 ! ! radiation absorbed by the layer-th ice layer 290 271 zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer) 291 272 END DO 292 273 END DO 293 274 294 ! Radiation transmitted below the ice 295 DO ji = kideb , kiut 296 fstbif_1d(ji) = fstbif_1d(ji) + & 297 zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) 275 DO ji = kideb, kiut ! Radiation transmitted below the ice 276 fstbif_1d(ji) = fstbif_1d(ji) + zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) 298 277 END DO 299 278 300 279 ! +++++ 301 280 ! just to check energy conservation 302 DO ji = kideb , kiut 303 zji = MOD( npb(ji) - 1, jpi ) + 1 304 zjj = ( npb(ji) - 1 ) / jpi + 1 305 fstroc(zji,zjj,jl) = & 306 zradtr_i(ji,nlay_i) 281 DO ji = kideb, kiut 282 ii = MOD( npb(ji) - 1, jpi ) + 1 283 ij = ( npb(ji) - 1 ) / jpi + 1 284 fstroc(ii,ij,jl) = zradtr_i(ji,nlay_i) 307 285 END DO 308 286 ! +++++ 309 287 310 288 DO layer = 1, nlay_i 311 DO ji = kideb 289 DO ji = kideb, kiut 312 290 radab(ji,layer) = zradab_i(ji,layer) 313 291 END DO … … 320 298 !------------------------------------------------------------------------------| 321 299 ! 322 ! Old surface temperature 323 DO ji = kideb, kiut 324 ztsuold(ji) = t_su_b(ji) ! temperature at the beg of iter pr. 325 ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter 326 t_su_b(ji) = MIN(t_su_b(ji),ztfs(ji)-0.00001) !necessary 327 zerrit(ji) = 1000.0 ! initial value of error 328 END DO 329 !RB Min global ?? 330 331 ! Old snow temperature 332 DO layer = 1, nlay_s 333 DO ji = kideb , kiut 334 ztsold(ji,layer) = t_s_b(ji,layer) 335 END DO 336 END DO 337 338 ! Old ice temperature 339 DO layer = 1, nlay_i 340 DO ji = kideb , kiut 341 ztiold(ji,layer) = t_i_b(ji,layer) 342 END DO 343 END DO 344 345 nconv = 0 ! number of iterations 346 zerritmax = 1000.0 ! maximal value of error on all points 347 348 DO WHILE ((zerritmax > maxer_i_thd).AND.(nconv < nconv_i_thd)) 349 350 nconv = nconv+1 351 300 DO ji = kideb, kiut ! Old surface temperature 301 ztsuold (ji) = t_su_b(ji) ! temperature at the beg of iter pr. 302 ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter 303 t_su_b (ji) = MIN( t_su_b(ji), ztfs(ji)-0.00001 ) ! necessary 304 zerrit (ji) = 1000._wp ! initial value of error 305 END DO 306 307 DO layer = 1, nlay_s ! Old snow temperature 308 DO ji = kideb , kiut 309 ztsold(ji,layer) = t_s_b(ji,layer) 310 END DO 311 END DO 312 313 DO layer = 1, nlay_i ! Old ice temperature 314 DO ji = kideb , kiut 315 ztiold(ji,layer) = t_i_b(ji,layer) 316 END DO 317 END DO 318 319 nconv = 0 ! number of iterations 320 zerritmax = 1000._wp ! maximal value of error on all points 321 322 DO WHILE ( zerritmax > maxer_i_thd .AND. nconv < nconv_i_thd ) 323 ! 324 nconv = nconv + 1 352 325 ! 353 326 !------------------------------------------------------------------------------| … … 355 328 !------------------------------------------------------------------------------| 356 329 ! 357 IF ( thcon_i_swi .EQ. 0 ) THEN 358 ! Untersteiner (1964) formula 330 IF( thcon_i_swi == 0 ) THEN ! Untersteiner (1964) formula 359 331 DO ji = kideb , kiut 360 332 ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / & … … 362 334 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 363 335 END DO 364 ENDIF365 366 IF ( thcon_i_swi .EQ. 1 ) THEN367 ! Pringle et al formula included,368 ! 2.11 + 0.09 S/T - 0.011.T369 DO ji = kideb , kiut370 ztcond_i(ji,0) = rcdic + 0.09*s_i_b(ji,1) / &371 MIN(-zeps,t_i_b(ji,1)-rtt) - &372 0.011* ( t_i_b(ji,1) - rtt )373 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin)374 END DO375 ENDIF376 377 IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner378 336 DO layer = 1, nlay_i-1 379 337 DO ji = kideb , kiut … … 406 364 ENDIF 407 365 408 IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle 409 DO ji = kideb , kiut 410 ztcond_i(ji,nlay_i) = rcdic + 0.09*s_i_b(ji,nlay_i) / & 411 MIN(-zeps,t_bo_b(ji)-rtt) - & 412 0.011* ( t_bo_b(ji) - rtt ) 413 ztcond_i(ji,nlay_i) = MAX(ztcond_i(ji,nlay_i),zkimin) 366 IF( thcon_i_swi == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 367 DO ji = kideb , kiut 368 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -zeps, t_i_b(ji,1)-rtt ) & 369 & - 0.011_wp * ( t_i_b(ji,1) - rtt ) 370 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 371 END DO 372 DO layer = 1, nlay_i-1 373 DO ji = kideb , kiut 374 ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) & 375 & / MIN(-2.0*zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1)-2.0*rtt) & 376 & - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt ) 377 ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 378 END DO 379 END DO 380 DO ji = kideb , kiut 381 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-zeps,t_bo_b(ji)-rtt) & 382 & - 0.011_wp * ( t_bo_b(ji) - rtt ) 383 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 414 384 END DO 415 385 ENDIF … … 732 702 733 703 ! surface temperature 734 isnow(ji) 735 ztsuoldit(ji) 704 isnow(ji) = INT(1.0-max(0.0,sign(1.0,-ht_s_b(ji)))) 705 ztsuoldit(ji) = t_su_b(ji) 736 706 IF (t_su_b(ji) .LT. ztfs(ji)) & 737 t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* & 738 ( isnow(ji)*t_s_b(ji,1) + & 739 (1.0-isnow(ji))*t_i_b(ji,1) ) ) / & 740 zdiagbis(ji,numeqmin(ji)) 707 t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( isnow(ji)*t_s_b(ji,1) & 708 & + (1.0-isnow(ji))*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 741 709 END DO 742 710 ! … … 748 716 ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 749 717 DO ji = kideb , kiut 750 t_su_b(ji) = MAX(MIN(t_su_b(ji),ztfs(ji)),190.0)751 zerrit(ji) = ABS(t_su_b(ji)-ztsuoldit(ji))718 t_su_b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp ) 719 zerrit(ji) = ABS( t_su_b(ji) - ztsuoldit(ji) ) 752 720 END DO 753 721 754 722 DO layer = 1, nlay_s 755 723 DO ji = kideb , kiut 756 zji = MOD( npb(ji) - 1, jpi ) + 1 757 zjj = ( npb(ji) - 1 ) / jpi + 1 758 t_s_b(ji,layer) = MAX(MIN(t_s_b(ji,layer),rtt),190.0) 759 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) & 760 - ztstemp(ji,layer))) 724 ii = MOD( npb(ji) - 1, jpi ) + 1 725 ij = ( npb(ji) - 1 ) / jpi + 1 726 t_s_b(ji,layer) = MAX( MIN( t_s_b(ji,layer), rtt ), 190._wp ) 727 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer))) 761 728 END DO 762 729 END DO … … 764 731 DO layer = 1, nlay_i 765 732 DO ji = kideb , kiut 766 ztmelt_i = -tmut*s_i_b(ji,layer) +rtt767 t_i_b(ji,layer) 768 zerrit(ji) 733 ztmelt_i = -tmut * s_i_b(ji,layer) + rtt 734 t_i_b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i),190.0) 735 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer))) 769 736 END DO 770 737 END DO 771 738 772 739 ! Compute spatial maximum over all errors 773 ! note that this could be optimized substantially by iterating only 774 ! the non-converging points 775 zerritmax = 0.0 776 DO ji = kideb , kiut 777 zerritmax = MAX(zerritmax,zerrit(ji)) 778 END DO 779 IF( lk_mpp ) CALL mpp_max(zerritmax, kcom=ncomm_ice) 740 ! note that this could be optimized substantially by iterating only the non-converging points 741 zerritmax = 0._wp 742 DO ji = kideb, kiut 743 zerritmax = MAX( zerritmax, zerrit(ji) ) 744 END DO 745 IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice ) 780 746 781 747 END DO ! End of the do while iterative procedure … … 787 753 788 754 ! 789 !-------------------------------------------------------------------------- 790 ! 11) Fluxes at the interfaces | 791 !-------------------------------------------------------------------------- 792 ! 755 !-------------------------------------------------------------------------! 756 ! 11) Fluxes at the interfaces ! 757 !-------------------------------------------------------------------------! 793 758 DO ji = kideb, kiut 794 ! update of latent heat fluxes 795 qla_ice_1d (ji) = qla_ice_1d (ji) + & 796 dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) 797 798 ! surface ice conduction flux 799 isnow(ji) = int(1.0-max(0.0,sign(1.0,-ht_s_b(ji)))) 800 fc_su(ji) = - isnow(ji)*zkappa_s(ji,0)*zg1s*(t_s_b(ji,1) - & 801 t_su_b(ji)) & 802 - (1.0-isnow(ji))*zkappa_i(ji,0)*zg1* & 803 (t_i_b(ji,1) - t_su_b(ji)) 804 805 ! bottom ice conduction flux 806 fc_bo_i(ji) = - zkappa_i(ji,nlay_i)* & 807 ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 808 759 ! ! update of latent heat fluxes 760 qla_ice_1d (ji) = qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) 761 ! ! surface ice conduction flux 762 isnow(ji) = INT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) ) 763 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) & 764 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji)) 765 ! ! bottom ice conduction flux 766 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 809 767 END DO 810 768 … … 812 770 ! Heat conservation ! 813 771 !-------------------------! 814 IF ( con_i ) THEN 815 772 IF( con_i ) THEN 816 773 DO ji = kideb, kiut 817 774 ! Upper snow value 818 fc_s(ji,0) = - isnow(ji)* & 819 zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - & 820 t_su_b(ji) ) 775 fc_s(ji,0) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) ) 821 776 ! Bott. snow value 822 fc_s(ji,1) = - isnow(ji)* & 823 zkappa_s(ji,1) * ( t_i_b(ji,1) - & 824 t_s_b(ji,1) ) 825 END DO 826 827 ! Upper ice layer 828 DO ji = kideb, kiut 777 fc_s(ji,1) = - isnow(ji)* zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) ) 778 END DO 779 DO ji = kideb, kiut ! Upper ice layer 829 780 fc_i(ji,0) = - isnow(ji) * & ! interface flux if there is snow 830 781 ( zkappa_i(ji,0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) & … … 832 783 zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not 833 784 END DO 834 835 ! Internal ice layers 836 DO layer = 1, nlay_i - 1 785 DO layer = 1, nlay_i - 1 ! Internal ice layers 837 786 DO ji = kideb, kiut 838 fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - & 839 t_i_b(ji,layer) ) 840 zji = MOD( npb(ji) - 1, jpi ) + 1 841 zjj = ( npb(ji) - 1 ) / jpi + 1 842 END DO 843 END DO 844 845 ! Bottom ice layers 846 DO ji = kideb, kiut 847 fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i)* & 848 ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 849 END DO 850 787 fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) 788 ii = MOD( npb(ji) - 1, jpi ) + 1 789 ij = ( npb(ji) - 1 ) / jpi + 1 790 END DO 791 END DO 792 DO ji = kideb, kiut ! Bottom ice layers 793 fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 794 END DO 851 795 ENDIF 852 796 ! 853 797 END SUBROUTINE lim_thd_dif 854 798 855 799 #else 856 !!====================================================================== 857 !! *** MODULE limthd_dif *** 858 !! no sea ice model 859 !!====================================================================== 800 !!---------------------------------------------------------------------- 801 !! Dummy Module No LIM-3 sea-ice model 802 !!---------------------------------------------------------------------- 860 803 CONTAINS 861 804 SUBROUTINE lim_thd_dif ! Empty routine
Note: See TracChangeset
for help on using the changeset viewer.
