Changeset 8522
- Timestamp:
- 2017-09-14T17:52:02+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO
- Files:
-
- 7 edited
-
LIM_SRC_3/icedia.F90 (modified) (1 diff)
-
LIM_SRC_3/icethd.F90 (modified) (2 diffs)
-
LIM_SRC_3/icethd_dh.F90 (modified) (3 diffs)
-
LIM_SRC_3/icethd_dif.F90 (modified) (36 diffs)
-
LIM_SRC_3/icevar.F90 (modified) (4 diffs)
-
LIM_SRC_3/icewri.F90 (modified) (3 diffs)
-
OPA_SRC/DOM/phycst.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90
r8518 r8522 19 19 USE phycst ! physical constant 20 20 USE daymod ! model calendar 21 USE sbc_oce , ONLY : sfx ! surface boundary condition: ocean fields21 USE sbc_oce , ONLY : sfx, nn_fsbc ! surface boundary condition: ocean fields 22 22 USE icerst ! ice restart 23 23 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
r8518 r8522 280 280 !!------------------------------------------------------------------- 281 281 INTEGER :: ji, jk ! dummy loop indices 282 REAL(wp) :: ztmelts, z 1_2cp, zbbb, zccc ! local scalar282 REAL(wp) :: ztmelts, zbbb, zccc ! local scalar 283 283 !!------------------------------------------------------------------- 284 284 ! Recover ice temperature 285 z1_2cp = 1._wp / ( 2._wp * cpic )286 285 DO jk = 1, nlay_i 287 286 DO ji = 1, nidx … … 290 289 zbbb = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus 291 290 zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) ) 292 t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * z1_2cp291 t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpic 293 292 294 293 ! mask temperature -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
r8514 r8522 78 78 REAL(wp) :: zgrr ! bottom growth rate 79 79 REAL(wp) :: zt_i_new ! bottom formation temperature 80 REAL(wp) :: z1_rho ! 1/(rhosn+rau0-rhoic) 80 81 81 82 REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean … … 596 597 ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, 597 598 ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 599 z1_rho = 1._wp / ( rhosn+rau0-rhoic ) 598 600 DO ji = 1, nidx 599 601 ! 600 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ))602 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) * z1_rho ) 601 603 602 604 ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) … … 646 648 e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 647 649 ! recalculate t_s_1d from e_s_1d 648 t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) / ( rhosn * cpic ) + lfus /cpic )650 t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhosn * r1_cpic + lfus * r1_cpic ) 649 651 END DO 650 652 END DO -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90
r8518 r8522 72 72 !! ice salinities : s_i_1d 73 73 !! number of layers in the ice/snow: nlay_i, nlay_s 74 !! profile of the ice/snow layers : z_i, z_s75 74 !! total ice/snow thickness : ht_i_1d, ht_s_1d 76 75 !!------------------------------------------------------------------ 77 INTEGER :: ji ! spatial loop index 78 INTEGER :: ii, ij ! temporary dummy loop index 76 INTEGER :: ji, jk ! spatial loop index 79 77 INTEGER :: numeq ! current reference number of equation 80 INTEGER :: jk ! vertical dummy loop index81 78 INTEGER :: minnumeqmin, maxnumeqmax 82 79 INTEGER :: iconv ! number of iterations in iterative procedure … … 93 90 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 94 91 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C 92 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 95 93 REAL(wp) :: ztmelt_i ! ice melting temperature 96 94 REAL(wp) :: z1_hsu 97 95 REAL(wp) :: zdti_max ! current maximal error on temperature 98 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature99 96 REAL(wp) :: zcpi ! Ice specific heat 100 REAL(wp) :: zhi !101 97 REAL(wp) :: zhfx_err, zdq ! diag errors on heat 98 REAL(wp) :: zfac ! dummy factor 102 99 103 100 REAL(wp), DIMENSION(jpij) :: isnow ! switch for presence (1) or absence (0) of snow 104 101 REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration 105 REAL(wp), DIMENSION(jpij) :: zh_i 106 REAL(wp), DIMENSION(jpij) :: zh_s 102 REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness 103 REAL(wp), DIMENSION(jpij) :: zh_s, z1_h_s ! snow layer thickness 107 104 REAL(wp), DIMENSION(jpij) :: zfsw ! solar radiation absorbed at the surface 108 105 REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface 109 106 REAL(wp), DIMENSION(jpij) :: zf ! surface flux function 110 107 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 111 REAL(wp), DIMENSION(jpij) :: zdti ! current error on temperature112 108 REAL(wp), DIMENSION(jpij) :: zftrice ! solar radiation transmitted through the ice 113 109 114 REAL(wp), DIMENSION(jpij,nlay_i) :: z_i ! Vertical cotes of the layers in the ice 115 REAL(wp), DIMENSION(jpij,nlay_s) :: z_s ! Vertical cotes of the layers in the snow 116 REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Old temperature in the ice 117 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow 118 REAL(wp), DIMENSION(jpij,nlay_i) :: ztitemp ! Temporary temperature in the ice to check the convergence 119 REAL(wp), DIMENSION(jpij,nlay_s) :: ztstemp ! Temporary temperature in the snow to check the convergence 110 REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice 111 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow 112 REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence 113 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence 120 114 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity 121 115 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice … … 136 130 ! Mono-category 137 131 REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done 138 REAL(wp) :: zratio_s ! dummy factor139 REAL(wp) :: zratio_i ! dummy factor140 REAL(wp) :: zh_thres ! thickness thres. for G(h) computation141 132 REAL(wp) :: zhe ! dummy factor 142 REAL(wp) :: zkimean ! mean sea ice thermal conductivity 143 REAL(wp) :: zfac ! dummy factor 144 REAL(wp) :: zihe ! dummy factor 145 REAL(wp) :: zheshth ! dummy factor 133 REAL(wp) :: zcnd_i ! mean sea ice thermal conductivity 146 134 !!------------------------------------------------------------------ 147 135 … … 155 143 ! 1) Initialization ! 156 144 !------------------------------------------------------------------------------! 157 DO ji = 1 145 DO ji = 1, nidx 158 146 isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ! is there snow or not 159 147 ! layer thickness … … 161 149 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 162 150 END DO 163 164 !-------------------- 165 ! Ice / snow layers 166 !-------------------- 167 z_s(1:nidx,1) = zh_s(1:nidx) 168 z_i(1:nidx,1) = zh_i(1:nidx) 169 DO jk = 2, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 170 DO ji = 1 , nidx 171 z_s(ji,jk) = z_s(ji,jk-1) + zh_s(ji) 172 END DO 173 END DO 174 175 DO jk = 2, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 176 DO ji = 1 , nidx 177 z_i(ji,jk) = z_i(ji,jk-1) + zh_i(ji) 178 END DO 179 END DO 151 ! 152 WHERE( zh_i(1:nidx) >= epsi10 ) ; z1_h_i(1:nidx) = 1._wp / zh_i(1:nidx) 153 ELSEWHERE ; z1_h_i(1:nidx) = 0._wp 154 END WHERE 155 156 WHERE( zh_s(1:nidx) >= epsi10 ) ; z1_h_s(1:nidx) = 1._wp / zh_s(1:nidx) 157 ELSEWHERE ; z1_h_s(1:nidx) = 0._wp 158 END WHERE 159 ! 160 ! temperatures 161 ztsub (1:nidx) = t_su_1d(1:nidx) ! temperature at the previous iteration 162 ztsold (1:nidx,:) = t_s_1d(1:nidx,:) ! Old snow temperature 163 ztiold (1:nidx,:) = t_i_1d(1:nidx,:) ! Old ice temperature 164 t_su_1d(1:nidx) = MIN( t_su_1d(1:nidx), rt0 - ztsu_err ) ! necessary 180 165 ! 181 166 !------------------------------------------------------------------------------| 182 ! 2) Radiation 167 ! 2) Radiation | 183 168 !------------------------------------------------------------------------------| 184 169 ! 185 170 z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 186 DO ji = 1 171 DO ji = 1, nidx 187 172 !------------------- 188 173 ! Computation of i0 … … 196 181 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 197 182 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 198 z hi= MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )199 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + z hi* fr2_i0_1d(ji) )183 zfac = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) ) 184 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) ) 200 185 201 186 !------------------------------------------------------- … … 216 201 DO ji = 1, nidx 217 202 ! ! radiation transmitted below the layer-th snow layer 218 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * z _s(ji,jk) )203 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) ) 219 204 ! ! radiation absorbed by the layer-th snow layer 220 205 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) … … 226 211 DO ji = 1, nidx 227 212 ! ! radiation transmitted below the layer-th ice layer 228 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * z _i(ji,jk) )213 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) 229 214 ! ! radiation absorbed by the layer-th ice layer 230 215 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) … … 238 223 !------------------------------------------------------------------------------| 239 224 ! 240 ztsub (1:nidx) = t_su_1d(1:nidx) ! temperature at the previous iter 241 t_su_1d(1:nidx) = MIN( t_su_1d(1:nidx), rt0 - ztsu_err ) ! necessary 242 zdti (1:nidx) = 1000._wp ! initial value of error 243 ztsb (1:nidx,:) = t_s_1d(1:nidx,:) ! Old snow temperature 244 ztib (1:nidx,:) = t_i_1d(1:nidx,:) ! Old ice temperature 245 246 iconv = 0 ! number of iterations 247 zdti_max = 1000._wp ! maximal value of error on all points 248 225 iconv = 0 ! number of iterations 226 zdti_max = 1000._wp ! maximal value of error on all points 249 227 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) 250 228 ! 251 229 iconv = iconv + 1 252 230 ! 231 ztib(1:nidx,:) = t_i_1d(1:nidx,:) 232 ztsb(1:nidx,:) = t_s_1d(1:nidx,:) 233 ! 253 234 !------------------------------------------------------------------------------| 254 235 ! 4) Sea ice thermal conductivity | … … 256 237 ! 257 238 IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T 258 DO ji = 1 , nidx 259 ztcond_i(ji,0) = MAX( zkimin, rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) ) 239 ! 240 DO ji = 1, nidx 241 ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 242 ztcond_i(ji,nlay_i) = rcdic + zbeta * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) 260 243 END DO 261 244 DO jk = 1, nlay_i-1 262 DO ji = 1 263 ztcond_i(ji,jk) = MAX( zkimin, rcdic + zbeta* ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / &264 & MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0))245 DO ji = 1, nidx 246 ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 247 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 265 248 END DO 266 249 END DO 267 250 ! 268 251 ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 269 DO ji = 1 , nidx 270 ztcond_i(ji,0) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 271 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) ) 252 ! 253 DO ji = 1, nidx 254 ztcond_i(ji,0) = rcdic + 0.09_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 255 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 256 ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 257 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 272 258 END DO 273 259 DO jk = 1, nlay_i-1 274 DO ji = 1 , nidx 275 ztcond_i(ji,jk) = MAX( zkimin, rcdic + & 276 & 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & 277 & / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 ) & 278 & - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) ) 260 DO ji = 1, nidx 261 ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 262 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) & 263 & - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 279 264 END DO 280 265 END DO 281 DO ji = 1 , nidx 282 ztcond_i(ji,nlay_i) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 283 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) ) 284 END DO 266 ! 285 267 ENDIF 286 268 ztcond_i(1:nidx,:) = MAX( zkimin, ztcond_i(1:nidx,:) ) 287 269 ! 288 270 !------------------------------------------------------------------------------| … … 294 276 ! Fichefet and Morales Maqueda, JGR 1997 295 277 296 zghe( :) = 1._wp297 278 zghe(1:nidx) = 1._wp 279 298 280 SELECT CASE ( nn_monocat ) 299 281 300 CASE ( 1,3)282 CASE ( 1 , 3 ) 301 283 302 284 zepsilon = 0.1_wp 303 zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp304 305 285 DO ji = 1, nidx 306 286 307 287 ! Mean sea ice thermal conductivity 308 z kimean = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp )288 zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) 309 289 310 290 ! Effective thickness he (zhe) 311 zfac = 1._wp / ( rn_cnd_s + zkimean ) 312 zratio_s = rn_cnd_s * zfac 313 zratio_i = zkimean * zfac 314 zhe = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 291 zhe = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) 315 292 316 293 ! G(he) 317 rswitch = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) ) ! =0 if zhe < zh_thres, if > 318 zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) 294 IF( zhe >= zepsilon * 0.5_wp * EXP(1._wp) ) THEN 295 zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) ) 296 ENDIF 319 297 320 ! Impose G(he) < 2.321 zghe(ji) = MIN( zghe(ji), 2._wp )322 323 298 END DO 324 299 325 300 END SELECT 326 327 301 ! 328 302 !------------------------------------------------------------------------------| … … 331 305 ! 332 306 !--- Snow 333 DO ji = 1, nidx 334 zfac = 1. / MAX( epsi10 , zh_s(ji) ) 335 zkappa_s(ji,0) = zghe(ji) * rn_cnd_s * zfac 336 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * zfac 337 END DO 338 339 DO jk = 1, nlay_s-1 340 DO ji = 1 , nidx 341 zkappa_s(ji,jk) = zghe(ji) * 2.0 * rn_cnd_s / MAX( epsi10, 2.0 * zh_s(ji) ) 342 END DO 307 DO jk = 0, nlay_s-1 308 DO ji = 1, nidx 309 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 310 END DO 311 END DO 312 DO ji = 1, nidx ! Snow-ice interface 313 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 314 IF( zfac > epsi10 ) THEN 315 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 316 ELSE 317 zkappa_s(ji,nlay_s) = 0._wp 318 ENDIF 343 319 END DO 344 320 345 321 !--- Ice 346 DO jk = 1, nlay_i-1 347 DO ji = 1 , nidx 348 zkappa_i(ji,jk) = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) ) 349 END DO 350 END DO 351 352 !--- Snow-ice interface 353 DO ji = 1 , nidx 354 zfac = 1./ MAX( epsi10 , zh_i(ji) ) 355 zkappa_i(ji,0) = zghe(ji) * ztcond_i(ji,0) * zfac 356 zkappa_i(ji,nlay_i) = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 357 zkappa_s(ji,nlay_s) = zghe(ji) * zghe(ji) * 2.0 * rn_cnd_s * ztcond_i(ji,0) / & 358 & MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cnd_s * zh_i(ji) ) ) 359 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 360 END DO 361 322 DO jk = 0, nlay_i 323 DO ji = 1, nidx 324 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 325 END DO 326 END DO 327 DO ji = 1, nidx ! Snow-ice interface 328 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 329 END DO 362 330 ! 363 331 !------------------------------------------------------------------------------| … … 366 334 ! 367 335 DO jk = 1, nlay_i 368 DO ji = 1 , nidx 369 ztitemp(ji,jk) = t_i_1d(ji,jk) 370 zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 371 zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zcpi * zh_i(ji), epsi10 ) 336 DO ji = 1, nidx 337 zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 338 zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) 372 339 END DO 373 340 END DO 374 341 375 342 DO jk = 1, nlay_s 376 DO ji = 1 , nidx 377 ztstemp(ji,jk) = t_s_1d(ji,jk) 378 zeta_s(ji,jk) = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 ) 379 END DO 380 END DO 381 343 DO ji = 1, nidx 344 zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 345 END DO 346 END DO 382 347 ! 383 348 !------------------------------------------------------------------------------| … … 386 351 ! 387 352 IF ( ln_dqns_i ) THEN 388 DO ji = 1 353 DO ji = 1, nidx 389 354 ! update of the non solar flux according to the update in T_su 390 355 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) … … 392 357 ENDIF 393 358 394 ! Update incoming flux 395 DO ji = 1 , nidx 396 ! update incoming flux 397 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 398 & + qns_ice_1d(ji) ! non solar total flux (LWup, LWdw, SH, LH) 399 END DO 400 359 DO ji = 1, nidx 360 zf(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 361 END DO 401 362 ! 402 363 !------------------------------------------------------------------------------| … … 412 373 413 374 DO numeq=1,nlay_i+3 414 DO ji = 1 375 DO ji = 1, nidx 415 376 ztrid(ji,numeq,1) = 0. 416 377 ztrid(ji,numeq,2) = 0. … … 423 384 424 385 DO numeq = nlay_s + 2, nlay_s + nlay_i 425 DO ji = 1 386 DO ji = 1, nidx 426 387 jk = numeq - nlay_s - 1 427 388 ztrid(ji,numeq,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 428 389 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 429 390 ztrid(ji,numeq,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 430 zindterm(ji,numeq) = zti b(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)391 zindterm(ji,numeq) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 431 392 END DO 432 393 ENDDO 433 394 434 395 numeq = nlay_s + nlay_i + 1 435 DO ji = 1 396 DO ji = 1, nidx 436 397 !!ice bottom term 437 398 ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) 438 399 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 439 400 ztrid(ji,numeq,3) = 0.0 440 zindterm(ji,numeq) = zti b(ji,nlay_i) + zeta_i(ji,nlay_i) * &401 zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 441 402 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 442 403 ENDDO 443 404 444 405 445 DO ji = 1 406 DO ji = 1, nidx 446 407 IF ( ht_s_1d(ji) > 0.0 ) THEN 447 408 ! … … 456 417 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 457 418 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 458 zindterm(ji,numeq) = zts b(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)419 zindterm(ji,numeq) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 459 420 END DO 460 421 461 422 !!case of only one layer in the ice (ice equation is altered) 462 IF ( nlay_i .eq.1 ) THEN423 IF ( nlay_i == 1 ) THEN 463 424 ztrid(ji,nlay_s+2,3) = 0.0 464 425 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) … … 483 444 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 484 445 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 485 zindterm(ji,2) = zts b(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)446 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 486 447 487 448 ELSE … … 498 459 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 499 460 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 500 zindterm(ji,2) = zts b(ji,1) + zeta_s(ji,1) * &461 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 501 462 & ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 502 463 ENDIF … … 526 487 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 527 488 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 528 zindterm(ji,numeqmin(ji)+1)= zti b(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)489 zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 529 490 530 491 !!case of only one layer in the ice (surface & ice equations are altered) … … 538 499 ztrid(ji,numeqmin(ji)+1,3) = 0.0 539 500 540 zindterm(ji,numeqmin(ji)+1) = zti b(ji,1) + zeta_i(ji,1) * &501 zindterm(ji,numeqmin(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * & 541 502 & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 542 503 ENDIF … … 556 517 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 557 518 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 558 zindterm(ji,numeqmin(ji)) = zti b(ji,1) + zeta_i(ji,1) * &519 zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1) * & 559 520 & ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 560 521 … … 564 525 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 565 526 ztrid(ji,numeqmin(ji),3) = 0.0 566 zindterm(ji,numeqmin(ji)) = zti b(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &527 zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) & 567 528 & + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 568 529 ENDIF … … 572 533 573 534 END DO 574 575 535 ! 576 536 !------------------------------------------------------------------------------| … … 578 538 !------------------------------------------------------------------------------| 579 539 ! 580 581 540 ! Solve the tridiagonal system with Gauss elimination method. 582 541 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, … … 586 545 minnumeqmin = nlay_i+5 587 546 588 DO ji = 1 547 DO ji = 1, nidx 589 548 zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji)) 590 549 zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2) … … 594 553 595 554 DO jk = minnumeqmin+1, maxnumeqmax 596 DO ji = 1 555 DO ji = 1, nidx 597 556 numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 598 557 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3) / zdiagbis(ji,numeq-1) … … 601 560 END DO 602 561 603 DO ji = 1 562 DO ji = 1, nidx 604 563 ! ice temperatures 605 t_i_1d(ji,nlay_i) 564 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji)) 606 565 END DO 607 566 608 567 DO numeq = nlay_i + nlay_s, nlay_s + 2, -1 609 DO ji = 1 568 DO ji = 1, nidx 610 569 jk = numeq - nlay_s - 1 611 570 t_i_1d(ji,jk) = ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq) … … 613 572 END DO 614 573 615 DO ji = 1 574 DO ji = 1, nidx 616 575 ! snow temperatures 617 IF (ht_s_1d(ji) > 0._wp) &618 t_s_1d(ji,nlay_s) 619 & / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) )620 576 IF( ht_s_1d(ji) > 0._wp ) THEN 577 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 578 & / zdiagbis(ji,nlay_s+1) 579 ENDIF 621 580 ! surface temperature 622 isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) )623 581 ztsub(ji) = t_su_1d(ji) 624 IF( t_su_1d(ji) < rt0 ) &582 IF( t_su_1d(ji) < rt0 ) THEN 625 583 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) * & 626 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 584 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 585 ENDIF 627 586 END DO 628 587 ! … … 632 591 ! 633 592 ! check that nowhere it has started to melt 634 ! zdti(ji) is a measure of error, it has to be under zdti_bnd 635 DO ji = 1 , nidx 636 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , 190._wp ) 637 zdti (ji) = ABS( t_su_1d(ji) - ztsub(ji) ) 593 ! zdti_max is a measure of error, it has to be under zdti_bnd 594 zdti_max = 0._wp 595 DO ji = 1, nidx 596 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 597 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 638 598 END DO 639 599 640 600 DO jk = 1, nlay_s 641 DO ji = 1 642 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), 190._wp)643 zdti (ji) = MAX( zdti(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )601 DO ji = 1, nidx 602 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 603 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 644 604 END DO 645 605 END DO 646 606 647 607 DO jk = 1, nlay_i 648 DO ji = 1 608 DO ji = 1, nidx 649 609 ztmelt_i = -tmut * s_i_1d(ji,jk) + rt0 650 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp )651 zdti (ji) = MAX( zdti(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) )610 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 611 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 652 612 END DO 653 613 END DO … … 655 615 ! Compute spatial maximum over all errors 656 616 ! note that this could be optimized substantially by iterating only the non-converging points 657 zdti_max = 0._wp658 DO ji = 1, nidx659 zdti_max = MAX( zdti_max, zdti(ji) )660 END DO661 617 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 662 618 663 619 END DO ! End of the do while iterative procedure 664 665 ! MV SIMIP 2016666 !--- Snow-ice interfacial temperature (diagnostic SIMIP)667 DO ji = 1, nidx668 zfac = 1. / MAX( epsi10 , rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) )669 IF( zh_s(ji) >= 1.e-3 ) THEN670 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + &671 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac672 ELSE673 t_si_1d(ji) = t_su_1d(ji)674 ENDIF675 END DO676 ! END MV SIMIP 2016677 620 678 621 IF( ln_icectl .AND. lwp ) THEN … … 687 630 DO ji = 1, nidx 688 631 ! ! surface ice conduction flux 689 isnow(ji) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )690 632 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 691 633 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) … … 693 635 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 694 636 END DO 695 696 ! MV SIMIP 2016697 !--- Conduction fluxes (positive downwards)698 DO ji = 1, nidx699 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)700 diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji) * a_i_1d(ji) / at_i_1d(ji)701 END DO702 ! END MV SIMIP 2016703 637 704 638 ! --- computes sea ice energy of melting compulsory for icethd_dh --- ! … … 730 664 731 665 END DO 666 667 ! --- Diagnostics SIMIP --- ! 668 DO ji = 1, nidx 669 !--- Conduction fluxes (positive downwards) 670 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 671 diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji) * a_i_1d(ji) / at_i_1d(ji) 672 673 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 674 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 675 IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 676 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 677 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 678 ELSE 679 t_si_1d(ji) = t_su_1d(ji) 680 ENDIF 681 END DO 732 682 ! 733 683 END SUBROUTINE ice_thd_dif … … 748 698 DO jk = 1, nlay_i ! Sea ice energy of melting 749 699 DO ji = 1, nidx 750 ztmelts = - tmut * s_i_1d(ji,jk) + rt0751 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point752 ! (sometimes dif scheme produces abnormally high temperatures)753 e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )&754 & + lfus * ( 1. 0 - ( ztmelts-rt0 )/ ( t_i_1d(ji,jk) - rt0 ) ) &755 & - rcp * ( ztmelts-rt0 ) )700 ztmelts = - tmut * s_i_1d(ji,jk) 701 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point 702 ! (sometimes dif scheme produces abnormally high temperatures) 703 e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) & 704 & + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) & 705 & - rcp * ztmelts ) 756 706 END DO 757 707 END DO -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
r8517 r8522 150 150 !!------------------------------------------------------------------ 151 151 INTEGER :: ji, jj, jk, jl ! dummy loop indices 152 REAL(wp) :: ze_i , z1_cp, z1_2cp! local scalars152 REAL(wp) :: ze_i ! local scalars 153 153 REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - - 154 154 REAL(wp) :: zhmax, z1_zhmax, zsm_i ! - - … … 193 193 194 194 !------------------- 195 ! Ice temperature [K] (with a minimum value (rt0 - 100.) imposed everywhere)195 ! Ice temperature [K] (with a minimum value (rt0 - 100.)) 196 196 !------------------- 197 197 zlay_i = REAL( nlay_i , wp ) ! number of layers 198 z1_2cp = 1._wp / ( 2._wp * cpic )199 198 DO jl = 1, jpl 200 199 DO jk = 1, nlay_i … … 208 207 zbbb = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus 209 208 zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) ) 210 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * z1_2cp, ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts209 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpic , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts 211 210 ! 212 211 ELSE !--- no ice 213 t_i(ji,jj,jk,jl) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C) 214 !!clem: I think we should impose rt0 instead 212 t_i(ji,jj,jk,jl) = rt0 215 213 ENDIF 216 214 END DO … … 220 218 221 219 !-------------------- 222 ! Snow temperature [K] (with a minimum value (rt0 - 100.) imposed everywhere)220 ! Snow temperature [K] (with a minimum value (rt0 - 100.)) 223 221 !-------------------- 224 222 zlay_s = REAL( nlay_s , wp ) 225 z1_cp = 1._wp / cpic226 223 DO jk = 1, nlay_s 227 224 WHERE( v_s(:,:,:) > epsi20 ) !--- icy area 228 t_s(:,:,jk,:) = MAX( -100._wp , MIN( z1_cp* ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0225 t_s(:,:,jk,:) = MAX( -100._wp , MIN( r1_cpic * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0 229 226 ELSEWHERE !--- no ice 230 t_s(:,:,jk,:) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C)227 t_s(:,:,jk,:) = rt0 231 228 END WHERE 232 229 END DO -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90
r8518 r8522 179 179 ENDIF 180 180 181 ztmp = rday /rhosn181 ztmp = rday * r1_rhosn 182 182 CALL iom_put( "vfxspr" , wfx_spr * ztmp ) ! precip (snow) 183 183 CALL iom_put( "vfxsnw" , wfx_snw * ztmp ) ! total snw growth/melt … … 237 237 ! Add-ons for SIMIP 238 238 !-------------------------------- 239 zrho1 = ( rau0 - rhoic ) / rau0; zrho2 = rhosn /rau0239 zrho1 = ( rau0 - rhoic ) * r1_rau0; zrho2 = rhosn * r1_rau0 240 240 241 241 IF ( iom_use( "icepres" ) ) CALL iom_put( "icepres" , zswi(:,:) ) ! Ice presence (1 or 0) … … 278 278 279 279 IF ( iom_use( "dmsspr" ) ) CALL iom_put( "dmsspr" , - wfx_spr ) ! Snow mass change through snow fall 280 IF ( iom_use( "dmsssi" ) ) CALL iom_put( "dmsssi" , wfx_sni*rhosn /rhoic) ! Snow mass change through snow-to-ice conversion280 IF ( iom_use( "dmsssi" ) ) CALL iom_put( "dmsssi" , wfx_sni*rhosn*r1_rhoic ) ! Snow mass change through snow-to-ice conversion 281 281 282 282 IF ( iom_use( "dmsmel" ) ) CALL iom_put( "dmsmel" , - wfx_snw_sum ) ! Snow mass change through melt -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
r8233 r8522 91 91 REAL(wp), PUBLIC :: r1_rhoic !: 1 / rhoic 92 92 REAL(wp), PUBLIC :: r1_rhosn !: 1 / rhosn 93 REAL(wp), PUBLIC :: r1_cpic !: 1 / cpic 93 94 #endif 94 95 !!---------------------------------------------------------------------- … … 159 160 r1_rhoic = 1._wp / rhoic 160 161 r1_rhosn = 1._wp / rhosn 162 r1_cpic = 1._wp / cpic 161 163 #endif 162 164 IF(lwp) THEN
Note: See TracChangeset
for help on using the changeset viewer.
