Changeset 2715 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
- Timestamp:
- 2011-03-30T17:58:35+02:00 (13 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r2528 r2715 7 7 !! ! 2005-06 (M. Vancoppenolle) 3D version 8 8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif & rdmicif 9 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 9 10 !!---------------------------------------------------------------------- 10 11 #if defined key_lim3 … … 17 18 USE phycst ! physical constants (OCE directory) 18 19 USE sbc_oce ! Surface boundary condition: ocean fields 19 USE ice 20 USE par_ice 21 USE thd_ice 22 USE in_out_manager 23 USE lib_mpp 20 USE ice ! LIM variables 21 USE par_ice ! LIM parameters 22 USE thd_ice ! LIM thermodynamics 23 USE wrk_nemo ! workspace manager 24 USE in_out_manager ! I/O manager 25 USE lib_mpp ! MPP library 24 26 25 27 IMPLICIT NONE … … 35 37 36 38 !!---------------------------------------------------------------------- 37 !! NEMO/LIM3 3.3, UCL - NEMO Consortium (2010)39 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 38 40 !! $Id$ 39 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)41 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 40 42 !!---------------------------------------------------------------------- 41 42 43 CONTAINS 43 44 44 SUBROUTINE lim_thd_dh( kideb,kiut,jl)45 SUBROUTINE lim_thd_dh( kideb, kiut, jl ) 45 46 !!------------------------------------------------------------------ 46 47 !! *** ROUTINE lim_thd_dh *** … … 75 76 INTEGER :: i_ice_switch ! ice thickness above a certain treshold or not 76 77 INTEGER :: iter 77 78 REAL(wp) :: zzfmass_i, zzfmass_s ! temporary scalar 79 REAL(wp) :: zhsnew, zihgnew, ztmelts ! temporary scalar 78 INTEGER :: num_iter_max, numce_dh 79 80 REAL(wp) :: meance_dh 81 REAL(wp) :: zzfmass_i, zihgnew ! local scalar 82 REAL(wp) :: zzfmass_s, zhsnew, ztmelts ! local scalar 80 83 REAL(wp) :: zhn, zdhcf, zdhbf, zhni, zhnfi, zihg ! 81 REAL(wp) :: zdhnm, zhnnew, zhisn, zihic 84 REAL(wp) :: zdhnm, zhnnew, zhisn, zihic, zzc ! 82 85 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 83 86 REAL(wp) :: zds ! increment of bottom ice salinity … … 89 92 REAL(wp) :: zgrr ! bottom growth rate 90 93 REAL(wp) :: ztform ! bottom formation temperature 91 92 REAL(wp), DIMENSION(jpij) :: zh_i ! ice layer thickness 93 REAL(wp), DIMENSION(jpij) :: zh_s ! snow layer thickness 94 REAL(wp), DIMENSION(jpij) :: ztfs ! melting point 95 REAL(wp), DIMENSION(jpij) :: zhsold ! old snow thickness 96 REAL(wp), DIMENSION(jpij) :: zqprec ! energy of fallen snow 97 REAL(wp), DIMENSION(jpij) :: zqfont_su ! incoming, remaining surface energy 98 REAL(wp), DIMENSION(jpij) :: zqfont_bo ! incoming, bottom energy 99 REAL(wp), DIMENSION(jpij) :: z_f_surf ! surface heat for ablation 100 REAL(wp), DIMENSION(jpij) :: zhgnew ! new ice thickness 101 REAL(wp), DIMENSION(jpij) :: zfmass_i ! 102 103 REAL(wp), DIMENSION(jpij) :: zdh_s_mel ! snow melt 104 REAL(wp), DIMENSION(jpij) :: zdh_s_pre ! snow precipitation 105 REAL(wp), DIMENSION(jpij) :: zdh_s_sub ! snow sublimation 106 REAL(wp), DIMENSION(jpij) :: zfsalt_melt ! salt flux due to ice melt 107 108 REAL(wp) , DIMENSION(jpij,jkmax) :: zdeltah 109 110 ! Pathological cases 111 REAL(wp), DIMENSION(jpij) :: zfdt_init ! total incoming heat for ice melt 112 REAL(wp), DIMENSION(jpij) :: zfdt_final ! total remaing heat for ice melt 113 REAL(wp), DIMENSION(jpij) :: zqt_i ! total ice heat content 114 REAL(wp), DIMENSION(jpij) :: zqt_s ! total snow heat content 115 REAL(wp), DIMENSION(jpij) :: zqt_dummy ! dummy heat content 116 94 ! 95 REAL(wp), POINTER, DIMENSION(:) :: zh_i, ztfs , zqfont_su, zqprec , zhgnew 96 REAL(wp), POINTER, DIMENSION(:) :: zh_s, zhsold, zqfont_bo, z_f_surf, zfmass_i 97 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel, zdh_s_sub , zfdt_init , zqt_i, zqt_dummy, zdq_i 98 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_pre, zfsalt_melt, zfdt_final, zqt_s, zfbase , zinnermelt 99 ! 100 REAL(wp), DIMENSION(jpij,jkmax) :: zdeltah 117 101 REAL(wp), DIMENSION(jpij,jkmax) :: zqt_i_lay ! total ice heat content 118 119 ! Heat conservation120 INTEGER :: num_iter_max, numce_dh121 REAL(wp) :: meance_dh122 INTEGER , DIMENSION(jpij) :: innermelt123 REAL(wp), DIMENSION(jpij) :: zfbase, zdq_i124 102 !!------------------------------------------------------------------ 125 103 126 zfsalt_melt(:) = 0.0 127 ftotal_fin(:) = 0.0 128 zfdt_init(:) = 0.0 129 zfdt_final(:) = 0.0 104 IF( wrk_in_use(1, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22) ) THEN 105 CALL ctl_stop('lim_thd_dh: requestead workspace arrays unavailable') ; RETURN 106 ENDIF 107 ! Set-up pointers to sub-arrays of workspace arrays 108 zh_i => wrk_1d_1 (1:jpij) ! ice layer thickness 109 zh_s => wrk_1d_2 (1:jpij) ! snow layer thickness 110 ztfs => wrk_1d_3 (1:jpij) ! melting point 111 zhsold => wrk_1d_4 (1:jpij) ! old snow thickness 112 zqprec => wrk_1d_5 (1:jpij) ! energy of fallen snow 113 zqfont_su => wrk_1d_6 (1:jpij) ! incoming, remaining surface energy 114 zqfont_bo => wrk_1d_7 (1:jpij) ! incoming, bottom energy 115 z_f_surf => wrk_1d_8 (1:jpij) ! surface heat for ablation 116 zhgnew => wrk_1d_9 (1:jpij) ! new ice thickness 117 zfmass_i => wrk_1d_10(1:jpij) ! 118 ! 119 zdh_s_mel => wrk_1d_11(1:jpij) ! snow melt 120 zdh_s_pre => wrk_1d_12(1:jpij) ! snow precipitation 121 zdh_s_sub => wrk_1d_13(1:jpij) ! snow sublimation 122 zfsalt_melt => wrk_1d_14(1:jpij) ! salt flux due to ice melt 123 ! 124 ! ! Pathological cases 125 zfdt_init => wrk_1d_15(1:jpij) ! total incoming heat for ice melt 126 zfdt_final => wrk_1d_16(1:jpij) ! total remaing heat for ice melt 127 zqt_i => wrk_1d_17(1:jpij) ! total ice heat content 128 zqt_s => wrk_1d_18(1:jpij) ! total snow heat content 129 zqt_dummy => wrk_1d_19(1:jpij) ! dummy heat content 130 131 zfbase => wrk_1d_20(1:jpij) 132 zdq_i => wrk_1d_21(1:jpij) 133 zinnermelt => wrk_1d_22(1:jpij) 134 135 zfsalt_melt(:) = 0._wp 136 ftotal_fin(:) = 0._wp 137 zfdt_init(:) = 0._wp 138 zfdt_final(:) = 0._wp 130 139 131 140 DO ji = kideb, kiut … … 138 147 !------------------------------------------------------------------------------! 139 148 ! 140 DO ji = kideb, kiut149 DO ji = kideb, kiut 141 150 isnow = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) ) 142 151 ztfs(ji) = isnow * rtt + ( 1.0 - isnow ) * rtt … … 146 155 END DO ! ji 147 156 148 zqfont_su (:) = 0.0149 zqfont_bo (:) = 0.0150 dsm_i_se_1d(:) = 0. 0151 dsm_i_si_1d(:) = 0. 0157 zqfont_su (:) = 0._wp 158 zqfont_bo (:) = 0._wp 159 dsm_i_se_1d(:) = 0._wp 160 dsm_i_si_1d(:) = 0._wp 152 161 ! 153 162 !------------------------------------------------------------------------------! … … 155 164 !------------------------------------------------------------------------------! 156 165 ! 157 ! Layer thickness 158 DO ji = kideb,kiut 166 DO ji = kideb, kiut ! Layer thickness 159 167 zh_i(ji) = ht_i_b(ji) / nlay_i 160 168 zh_s(ji) = ht_s_b(ji) / nlay_s 161 169 END DO 162 163 ! Total enthalpy of the snow 164 zqt_s(:) = 0.0 170 ! 171 zqt_s(:) = 0._wp ! Total enthalpy of the snow 165 172 DO jk = 1, nlay_s 166 DO ji = kideb, kiut173 DO ji = kideb, kiut 167 174 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 168 175 END DO 169 176 END DO 170 171 ! Total enthalpy of the ice 172 zqt_i(:) = 0.0 177 ! 178 zqt_i(:) = 0._wp ! Total enthalpy of the ice 173 179 DO jk = 1, nlay_i 174 DO ji = kideb,kiut 175 zqt_i(ji) = zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 176 zqt_i_lay(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 180 DO ji = kideb, kiut 181 zzc = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 182 zqt_i(ji) = zqt_i(ji) + zzc 183 zqt_i_lay(ji,jk) = zzc 177 184 END DO 178 185 END DO … … 201 208 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 202 209 END DO 203 zdh_s_mel(:) = 0. 0210 zdh_s_mel(:) = 0._wp 204 211 205 212 ! Melt of fallen snow … … 248 255 !-------------------------- 249 256 DO ji = kideb, kiut 250 dh_i_surf(ji) = 0. e0257 dh_i_surf(ji) = 0._wp 251 258 z_f_surf (ji) = zqfont_su(ji) / rdt_ice ! heat conservation test 252 zdq_i (ji) = 0. e0259 zdq_i (ji) = 0._wp 253 260 END DO ! ji 254 261 … … 267 274 ! 268 275 ! contribution to ice-ocean salt flux 269 zji = MOD( npb(ji) - 1 , jpi ) + 1270 zjj = ( npb(ji) - 1 ) / jpi + 1276 zji = MOD( npb(ji) - 1 , jpi ) + 1 277 zjj = ( npb(ji) - 1 ) / jpi + 1 271 278 zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji) ) * a_i_b(ji) & 272 279 & * MIN( zdeltah(ji,jk) , 0.e0 ) * rhoic / rdt_ice … … 278 285 ! !------------------- 279 286 numce_dh = 0 280 meance_dh = 0. e0287 meance_dh = 0._wp 281 288 DO ji = kideb, kiut 282 289 IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN … … 287 294 WRITE(numout,*) ' ALERTE heat loss for surface melt ' 288 295 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl 289 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)290 WRITE(numout,*) ' z_f_surf : ', z_f_surf(ji)291 WRITE(numout,*) ' zdq_i : ', zdq_i(ji)292 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)293 WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)294 WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)295 WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(ji)296 WRITE(numout,*) ' s_i_new : ', s_i_new(ji)297 WRITE(numout,*) ' sss_m : ', sss_m(zji,zjj)296 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 297 WRITE(numout,*) ' z_f_surf : ', z_f_surf(ji) 298 WRITE(numout,*) ' zdq_i : ', zdq_i(ji) 299 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 300 WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji) 301 WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji) 302 WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(ji) 303 WRITE(numout,*) ' s_i_new : ', s_i_new(ji) 304 WRITE(numout,*) ' sss_m : ', sss_m(zji,zjj) 298 305 ENDIF 299 306 END DO … … 440 447 ! 4.2 Basal melt 441 448 !---------------- 442 meance_dh = 0. 0449 meance_dh = 0._wp 443 450 numce_dh = 0 444 innermelt(:) = 0451 zinnermelt(:) = 0._wp 445 452 446 453 DO ji = kideb, kiut 447 454 ! heat convergence at the surface > 0 448 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0 ) THEN 449 455 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp ) THEN 450 456 s_i_new(ji) = s_i_b(ji,nlay_i) 451 457 zqfont_bo(ji) = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) 452 453 zfbase(ji) = zqfont_bo(ji) / rdt_ice ! heat conservation test 454 zdq_i(ji) = 0.e0 455 456 dh_i_bott(ji) = 0.e0 458 zfbase(ji) = zqfont_bo(ji) / rdt_ice ! heat conservation test 459 zdq_i(ji) = 0._wp 460 dh_i_bott(ji) = 0._wp 457 461 ENDIF 458 462 END DO … … 461 465 DO ji = kideb, kiut 462 466 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN 463 ztmelts 464 IF ( t_i_b(ji,jk) .GE.ztmelts ) THEN467 ztmelts = - tmut * s_i_b(ji,jk) + rtt 468 IF( t_i_b(ji,jk) >= ztmelts ) THEN 465 469 zdeltah(ji,jk) = - zh_i(ji) 466 470 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) 467 innermelt(ji) = 1471 zinnermelt(ji) = 1._wp 468 472 ELSE ! normal ablation 469 473 zdeltah(ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) … … 492 496 ENDIF 493 497 IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN 494 WRITE(numout,*) ' ALERTE heat loss for basal melt ' 495 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl 496 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 497 WRITE(numout,*) ' zfbase : ', zfbase(ji) 498 WRITE(numout,*) ' zdq_i : ', zdq_i(ji) 499 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 500 WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji) 501 WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji) 502 WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji) 503 WRITE(numout,*) ' s_i_new : ', s_i_new(ji) 504 WRITE(numout,*) ' sss_m : ', sss_m(zji,zjj) 498 WRITE(numout,*) ' ALERTE heat loss for basal melt : zji, zjj, jl :', zji, zjj, jl 499 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 500 WRITE(numout,*) ' zfbase : ', zfbase(ji) 501 WRITE(numout,*) ' zdq_i : ', zdq_i(ji) 502 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 503 WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji) 504 WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji) 505 WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(ji) 506 WRITE(numout,*) ' s_i_new : ', s_i_new(ji) 507 WRITE(numout,*) ' sss_m : ', sss_m(zji,zjj) 505 508 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 506 WRITE(numout,*) ' innermelt : ', innermelt(ji)509 WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) ) 507 510 ENDIF 508 511 ENDIF … … 687 690 688 691 ! Total ablation ! new lines added to debug 689 IF( ht_i_b(ji) <= 0. e0 ) a_i_b(ji) = 0.0692 IF( ht_i_b(ji) <= 0._wp ) a_i_b(ji) = 0._wp 690 693 691 694 ! diagnostic ( snow ice growth ) … … 695 698 ! 696 699 END DO !ji 697 700 ! 701 IF( wrk_not_released(1, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21) ) & 702 CALL ctl_stop('lim_thd_dh : failed to release workspace arrays') 703 ! 698 704 END SUBROUTINE lim_thd_dh 699 705
Note: See TracChangeset
for help on using the changeset viewer.
