Changeset 1218 for trunk/NEMO/LIM_SRC_2
- Timestamp:
- 2008-10-28T10:12:16+01:00 (16 years ago)
- Location:
- trunk/NEMO/LIM_SRC_2
- Files:
-
- 3 edited
-
limsbc_2.F90 (modified) (6 diffs)
-
limthd_2.F90 (modified) (24 diffs)
-
limthd_zdf_2.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_2/limsbc_2.F90
r1173 r1218 29 29 USE albedo ! albedo parameters 30 30 USE prtctl ! Print control 31 USE cpl_oasis3, ONLY : lk_cpl 31 32 32 33 IMPLICIT NONE … … 85 86 REAL(wp) :: zutau , zvtau ! lead fraction at U- & V-points 86 87 REAL(wp) :: zu_io , zv_io ! 2 components of the ice-ocean velocity 87 #if defined key_coupled 88 REAL(wp), DIMENSION(jpi,jpj) :: zalb ! albedo of ice under overcast sky 89 REAL(wp), DIMENSION(jpi,jpj) :: zalbp ! albedo of ice under clear sky 90 #endif 88 ! interface 2D --> 3D 89 REAL(wp), DIMENSION(jpi,jpj,1) :: zalb ! albedo of ice under overcast sky 90 REAL(wp), DIMENSION(jpi,jpj,1) :: zalbp ! albedo of ice under clear sky 91 REAL(wp), DIMENSION(jpi,jpj,1) :: zsist ! surface ice temperature (K) 92 REAL(wp), DIMENSION(jpi,jpj,1) :: zhicif ! ice thickness 93 REAL(wp), DIMENSION(jpi,jpj,1) :: zhsnif ! snow thickness 91 94 REAL(wp) :: zsang, zmod, zfm 92 95 REAL(wp), DIMENSION(jpi,jpj) :: ztio_u, ztio_v ! ocean stress below sea-ice … … 119 122 ifral = ( 1 - i1mfr * ( 1 - ial ) ) 120 123 ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv 124 125 !!$ zinda = 1.0 - AINT( pfrld(ji,jj) ) ! = 0. if pure ocean else 1. (at previous time) 126 !!$ 127 !!$ i1mfr = 1.0 - AINT( frld(ji,jj) ) ! = 0. if pure ocean else 1. (at current time) 128 !!$ 129 !!$ IF( phicif(ji,jj) <= 0. ) THEN ; ifvt = zinda ! = 1. if (snow and no ice at previous time) else 0. ??? 130 !!$ ELSE ; ifvt = 0. 131 !!$ ENDIF 132 !!$ 133 !!$ IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN ; idfr = 0. ! = 0. if lead fraction increases from previous to current 134 !!$ ELSE ; idfr = 1. 135 !!$ ENDIF 136 !!$ 137 !!$ iflt = zinda * (1 - i1mfr) * (1 - ifvt ) ! = 1. if ice (not only snow) at previous and pure ocean at current 138 !!$ 139 !!$ ial = ifvt * i1mfr + ( 1 - ifvt ) * idfr 140 !!$! snow no ice ice ice or nothing lead fraction increases 141 !!$! at previous now at previous 142 !!$! -> ice aera increases ??? -> ice aera decreases ??? 143 !!$ 144 !!$ iadv = ( 1 - i1mfr ) * zinda 145 !!$! pure ocean ice at 146 !!$! at current previous 147 !!$! -> = 1. if ice disapear between previous and current 148 !!$ 149 !!$ ifral = ( 1 - i1mfr * ( 1 - ial ) ) 150 !!$! ice at ??? 151 !!$! current 152 !!$! -> ??? 153 !!$ 154 !!$ ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv 155 !!$! ice disapear 156 !!$ 157 !!$ 158 121 159 ! computation the solar flux at ocean surface 122 zqsr = pfrld(ji,jj) * qsr(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) 160 #if defined key_coupled 161 zqsr = tqsr(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj) ) * ( 1.0 - pfrld(ji,jj) ) 162 #else 163 zqsr = pfrld(ji,jj) * qsr(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) 164 #endif 123 165 ! computation the non solar heat flux at ocean surface 124 166 zqns = - ( 1. - thcm(ji,jj) ) * zqsr & ! part of the solar energy used in leads … … 145 187 DO ji = 1, jpi 146 188 189 #if defined key_coupled 190 zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) ) & ! 191 & + rdmsnif(ji,jj) / rdt_ice ! freshwaterflux due to snow melting 192 #else 193 !!$ ! computing freshwater exchanges at the ice/ocean interface 194 !!$ zpme = - evap(ji,jj) * frld(ji,jj) & ! evaporation over oceanic fraction 195 !!$ & + tprecip(ji,jj) & ! total precipitation 196 !!$ & - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) ) & ! remov. snow precip over ice 197 !!$ & - rdmsnif(ji,jj) / rdt_ice ! freshwaterflux due to snow melting 147 198 ! computing freshwater exchanges at the ice/ocean interface 148 199 zemp = + emp(ji,jj) * frld(ji,jj) & ! e-p budget over open ocean fraction … … 151 202 & + rdmsnif(ji,jj) / rdt_ice ! freshwaterflux due to snow melting 152 203 ! ! ice-covered fraction: 204 #endif 153 205 154 206 ! computing salt exchanges at the ice/ocean interface … … 217 269 218 270 fr_i (:,:) = 1.0 - frld(:,:) ! sea-ice fraction 219 tn_ice(:,:) = sist(:,:) ! sea-ice surface temperature 220 221 #if defined key_coupled 222 !------------------------------------------------! 223 ! Computation of snow/ice and ocean albedo ! 224 !------------------------------------------------! 225 zalb (:,:) = 0.e0 226 zalbp (:,:) = 0.e0 227 228 CALL albedo_ice( sist, hicif, hsnif, zalbp, zalb ) 229 230 alb_ice(:,:) = 0.5 * zalbp(:,:) + 0.5 * zalb (:,:) ! Ice albedo (mean clear and overcast skys) 231 #endif 271 272 IF ( lk_cpl ) THEN 273 ! Ice surface temperature 274 tn_ice(:,:) = sist(:,:) ! sea-ice surface temperature 275 ! Computation of snow/ice and ocean albedo 276 ! INTERFACE 3D versus 2D 277 zsist (:,:,1) = sist (:,:) 278 zhicif(:,:,1) = hicif(:,:) ; zhsnif(:,:,1) = hsnif(:,:) 279 CALL albedo_ice( zsist, zhicif, zhsnif, zalbp, zalb ) 280 alb_ice(:,:) = 0.5 * ( zalbp(:,:,1) + zalb (:,:,1) ) ! Ice albedo (mean clear and overcast skys) 281 ENDIF 232 282 233 283 IF(ln_ctl) THEN -
trunk/NEMO/LIM_SRC_2/limthd_2.F90
r1156 r1218 7 7 !! 2.0 ! 02-07 (C. Ethe, G. Madec) F90 8 8 !! 2.0 ! 03-08 (C. Ethe) add lim_thd_init 9 !! - ! 08-2008 (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface 9 10 !!--------------------------------------------------------------------- 10 11 #if defined key_lim2 … … 15 16 !! lim_thd_init_2 : initialisation of sea-ice thermodynamic 16 17 !!---------------------------------------------------------------------- 17 !! * Modules used18 18 USE phycst ! physical constants 19 19 USE dom_oce ! ocean space and time domain variables … … 31 31 USE limtab_2 32 32 USE prtctl ! Print control 33 USE cpl_oasis3, ONLY : lk_cpl 33 34 34 35 IMPLICIT NONE … … 37 38 PUBLIC lim_thd_2 ! called by lim_step 38 39 39 REAL(wp) :: epsi20 = 1.e-20 , &! constant values40 & epsi16 = 1.e-16 , &41 & epsi04 = 1.e-04 , &42 & zzero = 0.e0 , &43 & zone = 1.e040 REAL(wp) :: epsi20 = 1.e-20 ! constant values 41 REAL(wp) :: epsi16 = 1.e-16 ! 42 REAL(wp) :: epsi04 = 1.e-04 ! 43 REAL(wp) :: rzero = 0.e0 ! 44 REAL(wp) :: rone = 1.e0 ! 44 45 45 46 !! * Substitutions … … 47 48 # include "vectopt_loop_substitute.h90" 48 49 !!-------- ------------------------------------------------------------- 49 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)50 !! NEMO/LIM 2.0, UCL-LOCEAN-IPSL (2008) 50 51 !! $Id$ 51 52 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 74 75 INTEGER, INTENT(in) :: kt ! number of iteration 75 76 !! 76 INTEGER :: ji, jj , &! dummy loop indices77 nbpb , &! nb of icy pts for thermo. cal.78 nbpac! nb of pts for lateral accretion77 INTEGER :: ji, jj ! dummy loop indices 78 INTEGER :: nbpb ! nb of icy pts for thermo. cal. 79 INTEGER :: nbpac ! nb of pts for lateral accretion 79 80 CHARACTER (len=22) :: charout 80 REAL(wp) :: & 81 zfric_umin = 5e-03 , & ! lower bound for the friction velocity 82 zfric_umax = 2e-02 ! upper bound for the friction velocity 83 REAL(wp) :: & 84 zinda , & ! switch for test. the val. of concen. 85 zindb, zindg , & ! switches for test. the val of arg 86 za , zh, zthsnice , & 87 zfric_u , & ! friction velocity 88 zfnsol , & ! total non solar heat 89 zfontn , & ! heat flux from snow thickness 90 zfntlat, zpareff ! test. the val. of lead heat budget 91 REAL(wp), DIMENSION(jpi,jpj) :: zhicifp, & ! ice thickness for outputs 92 & zqlbsbq ! link with lead energy budget qldif 81 REAL(wp) :: zfric_umin = 5e-03 ! lower bound for the friction velocity 82 REAL(wp) :: zfric_umax = 2e-02 ! upper bound for the friction velocity 83 REAL(wp) :: zinda ! switch for test. the val. of concen. 84 REAL(wp) :: zindb, zindg ! switches for test. the val of arg 85 REAL(wp) :: za , zh, zthsnice ! 86 REAL(wp) :: zfric_u ! friction velocity 87 REAL(wp) :: zfnsol ! total non solar heat 88 REAL(wp) :: zfontn ! heat flux from snow thickness 89 REAL(wp) :: zfntlat, zpareff ! test. the val. of lead heat budget 90 REAL(wp) :: zfi ! temporary scalar 91 REAL(wp), DIMENSION(jpi,jpj) :: zhicifp ! ice thickness for outputs 92 REAL(wp), DIMENSION(jpi,jpj) :: zqlbsbq ! link with lead energy budget qldif 93 93 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmsk ! working array 94 94 !!------------------------------------------------------------------- … … 100 100 !-------------------------------------------! 101 101 102 ! i est-ce utile? oui au moins en partie102 !!gm needed? yes at least for some of these arrays 103 103 rdvosif(:,:) = 0.e0 ! variation of ice volume at surface 104 104 rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom … … 114 114 zmsk (:,:,:) = 0.e0 115 115 116 DO jj = 1, jpj 117 DO ji = 1, jpi 118 hsnif(ji,jj) = hsnif(ji,jj) * MAX( zzero, SIGN( zone , hsnif(ji,jj) - epsi04 ) ) 119 END DO 120 END DO 121 122 IF(ln_ctl) CALL prt_ctl(tab2d_1=hsnif , clinfo1=' lim_thd: hsnif : ') 116 ! set to zero snow thickness smaller than epsi04 117 DO jj = 1, jpj 118 DO ji = 1, jpi 119 hsnif(ji,jj) = hsnif(ji,jj) * MAX( rzero, SIGN( rone , hsnif(ji,jj) - epsi04 ) ) 120 END DO 121 END DO 122 !!gm better coded (do not use SIGN...) 123 ! WHERE( hsnif(:,:) < epsi04 ) hsnif(:,:) = 0.e0 124 !!gm 125 126 IF(ln_ctl) CALL prt_ctl( tab2d_1=hsnif, clinfo1=' lim_thd: hsnif : ' ) 123 127 124 128 !-----------------------------------! … … 129 133 DO ji = 1, jpi 130 134 ! snow is transformed into ice if the original ice cover disappears. 131 zindg = tms(ji,jj) * MAX( zzero , SIGN( zone , -hicif(ji,jj) ) )135 zindg = tms(ji,jj) * MAX( rzero , SIGN( rone , -hicif(ji,jj) ) ) 132 136 hicif(ji,jj) = hicif(ji,jj) + zindg * rhosn * hsnif(ji,jj) / rau0 133 hsnif(ji,jj) = ( zone - zindg ) * hsnif(ji,jj) + zindg * hicif(ji,jj) * ( rau0 - rhoic ) / rhosn137 hsnif(ji,jj) = ( rone - zindg ) * hsnif(ji,jj) + zindg * hicif(ji,jj) * ( rau0 - rhoic ) / rhosn 134 138 dmgwi(ji,jj) = zindg * (1.0 - frld(ji,jj)) * rhoic * hicif(ji,jj) ! snow/ice mass 135 139 136 140 ! the lead fraction, frld, must be little than or equal to amax (ice ridging). 137 141 zthsnice = hsnif(ji,jj) + hicif(ji,jj) 138 zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) )139 za = zindb * MIN( zone, ( 1.0 - frld(ji,jj) ) * uscomi )142 zindb = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) ) 143 za = zindb * MIN( rone, ( 1.0 - frld(ji,jj) ) * uscomi ) 140 144 hsnif (ji,jj) = hsnif(ji,jj) * za 141 145 hicif (ji,jj) = hicif(ji,jj) * za 142 146 qstoif(ji,jj) = qstoif(ji,jj) * za 143 frld (ji,jj) = 1.0 - zindb * ( 1.0 - frld(ji,jj) ) / MAX( za 147 frld (ji,jj) = 1.0 - zindb * ( 1.0 - frld(ji,jj) ) / MAX( za, epsi20 ) 144 148 145 149 ! the in situ ice thickness, hicif, must be equal to or greater than hiclim. 146 zh = MAX( zone , zindb * hiclim / MAX( hicif(ji,jj), epsi20 ) )150 zh = MAX( rone , zindb * hiclim / MAX( hicif(ji,jj), epsi20 ) ) 147 151 hsnif (ji,jj) = hsnif(ji,jj) * zh 148 152 hicif (ji,jj) = hicif(ji,jj) * zh … … 153 157 154 158 IF(ln_ctl) THEN 155 CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif : ')156 CALL prt_ctl( tab2d_1=hsnif , clinfo1=' lim_thd: hsnif : ')157 CALL prt_ctl( tab2d_1=dmgwi , clinfo1=' lim_thd: dmgwi : ')158 CALL prt_ctl( tab2d_1=qstoif , clinfo1=' lim_thd: qstoif : ')159 CALL prt_ctl( tab2d_1=frld , clinfo1=' lim_thd: frld : ')159 CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif : ' ) 160 CALL prt_ctl( tab2d_1=hsnif , clinfo1=' lim_thd: hsnif : ' ) 161 CALL prt_ctl( tab2d_1=dmgwi , clinfo1=' lim_thd: dmgwi : ' ) 162 CALL prt_ctl( tab2d_1=qstoif, clinfo1=' lim_thd: qstoif : ' ) 163 CALL prt_ctl( tab2d_1=frld , clinfo1=' lim_thd: frld : ' ) 160 164 ENDIF 161 165 … … 175 179 DO ji = 1, jpi 176 180 zthsnice = hsnif(ji,jj) + hicif(ji,jj) 177 zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) )181 zindb = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) ) 178 182 pfrld(ji,jj) = frld(ji,jj) 179 zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )183 zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 180 184 181 185 ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 189 193 190 194 ! partial computation of the lead energy budget (qldif) 191 zfontn = ( sprecip(ji,jj) / rhosn ) * xlsn ! energy for melting 195 #if defined key_coupled 196 zfi = 1.0 - pfrld(ji,jj) 197 qldif(ji,jj) = tms(ji,jj) * rdt_ice & 198 & * ( ( qsr_tot(ji,jj) - qsr_ice(ji,jj) * zfi ) * ( 1.0 - thcm(ji,jj) ) & 199 & + ( qns_tot(ji,jj) - qns_ice(ji,jj) * zfi ) & 200 & + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) ) ) 201 #else 202 zfontn = ( sprecip(ji,jj) / rhosn ) * xlsn ! energy for melting solid precipitation 192 203 zfnsol = qns(ji,jj) ! total non solar flux over the ocean 193 204 qldif(ji,jj) = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) ) & 194 205 & + zfnsol + fdtcn(ji,jj) - zfontn & 195 206 & + ( 1.0 - zindb ) * fsbbq(ji,jj) ) & 196 & * frld(ji,jj) * rdt_ice 207 & * frld(ji,jj) * rdt_ice 208 !!$ qldif(ji,jj) = tms(ji,jj) * rdt_ice * frld(ji,jj) 209 !!$ & * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) ) & 210 !!$ & + qns(ji,jj) + fdtcn(ji,jj) - zfontn & 211 !!$ & + ( 1.0 - zindb ) * fsbbq(ji,jj) ) & 212 #endif 197 213 ! parlat : percentage of energy used for lateral ablation (0.0) 198 zfntlat = 1.0 - MAX( zzero , SIGN( zone , - qldif(ji,jj) ) )214 zfntlat = 1.0 - MAX( rzero , SIGN( rone , - qldif(ji,jj) ) ) 199 215 zpareff = 1.0 + ( parlat - 1.0 ) * zinda * zfntlat 200 216 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( (1.0 - frld(ji,jj)) * rdt_ice , epsi16 ) … … 243 259 !------------------------------------------------------------------------------------ 244 260 245 IF ( nbpb > 0) THEN246 261 IF( nbpb > 0 ) THEN 262 ! 247 263 ! put the variable in a 1-D array for thermodynamics process 248 264 CALL tab_2d_1d_2( nbpb, frld_1d (1:nbpb) , frld , jpi, jpj, npb(1:nbpb) ) … … 257 273 CALL tab_2d_1d_2( nbpb, fr2_i0_1d (1:nbpb) , fr2_i0 , jpi, jpj, npb(1:nbpb) ) 258 274 CALL tab_2d_1d_2( nbpb, qns_ice_1d (1:nbpb) , qns_ice , jpi, jpj, npb(1:nbpb) ) 259 #if ! defined key_coupled 260 CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb) , qla_ice , jpi, jpj, npb(1:nbpb) )261 CALL tab_2d_1d_2( nbpb, dqla_ice_1d(1:nbpb) , dqla_ice , jpi, jpj, npb(1:nbpb) )262 #endif 275 IF( .NOT. lk_cpl ) THEN 276 CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb) , qla_ice , jpi, jpj, npb(1:nbpb) ) 277 CALL tab_2d_1d_2( nbpb, dqla_ice_1d(1:nbpb) , dqla_ice , jpi, jpj, npb(1:nbpb) ) 278 ENDIF 263 279 CALL tab_2d_1d_2( nbpb, dqns_ice_1d(1:nbpb) , dqns_ice , jpi, jpj, npb(1:nbpb) ) 264 280 CALL tab_2d_1d_2( nbpb, tfu_1d (1:nbpb) , tfu , jpi, jpj, npb(1:nbpb) ) … … 271 287 CALL tab_2d_1d_2( nbpb, dmgwi_1d (1:nbpb) , dmgwi , jpi, jpj, npb(1:nbpb) ) 272 288 CALL tab_2d_1d_2( nbpb, qlbbq_1d (1:nbpb) , zqlbsbq , jpi, jpj, npb(1:nbpb) ) 273 289 ! 274 290 CALL lim_thd_zdf_2( 1, nbpb ) ! compute ice growth 275 291 ! 276 292 ! back to the geographic grid. 277 293 CALL tab_1d_2d_2( nbpb, frld , npb, frld_1d (1:nbpb) , jpi, jpj ) … … 295 311 CALL tab_1d_2d_2( nbpb, fdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj ) 296 312 CALL tab_1d_2d_2( nbpb, rdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj ) 297 298 299 ENDIF 300 301 302 ! Up-date sea ice thickness. 303 !--------------------------------- 313 ! 314 ENDIF 315 316 317 ! Up-date sea ice thickness 318 !-------------------------- 304 319 DO jj = 1, jpj 305 320 DO ji = 1, jpi 306 321 phicif(ji,jj) = hicif(ji,jj) 307 hicif(ji,jj) = hicif(ji,jj) * ( zone - MAX( zzero, SIGN( zone, - ( 1.0 - frld(ji,jj) ) ) ) )308 END DO 309 END DO 310 311 312 ! Tricky trick : add 2 to frld in the Southern Hemisphere.313 !-------------------------------------------------------- --322 hicif(ji,jj) = hicif(ji,jj) * ( rone - MAX( rzero, SIGN( rone, - ( 1.0 - frld(ji,jj) ) ) ) ) 323 END DO 324 END DO 325 326 327 ! Tricky trick : add 2 to frld in the Southern Hemisphere 328 !-------------------------------------------------------- 314 329 IF( fcor(1,1) < 0.e0 ) THEN 315 330 DO jj = 1, njeqm1 … … 321 336 322 337 323 ! 324 ! 338 ! Select points for lateral accretion (this occurs when heat exchange 339 ! between ice and ocean is negative; ocean losing heat) 325 340 !----------------------------------------------------------------- 326 341 nbpac = 0 … … 341 356 ENDIF 342 357 343 344 ! 345 ! If ocean gains heat do nothing ; otherwise, one performs lateral accretion 358 359 ! If ocean gains heat do nothing ; otherwise, one performs lateral accretion 346 360 !-------------------------------------------------------------------------------- 347 348 361 IF( nbpac > 0 ) THEN 349 362 ! 350 363 !...Put the variable in a 1-D array for lateral accretion 351 364 CALL tab_2d_1d_2( nbpac, frld_1d (1:nbpac) , frld , jpi, jpj, npac(1:nbpac) ) … … 361 374 CALL tab_2d_1d_2( nbpac, dvlbq_1d (1:nbpac) , fdvolif , jpi, jpj, npac(1:nbpac) ) 362 375 CALL tab_2d_1d_2( nbpac, tfu_1d (1:nbpac) , tfu , jpi, jpj, npac(1:nbpac) ) 363 364 ! call lateral accretion routine. 365 CALL lim_thd_lac_2( 1 , nbpac ) 366 376 ! 377 CALL lim_thd_lac_2( 1 , nbpac ) ! lateral accretion routine. 378 ! 367 379 ! back to the geographic grid 368 380 CALL tab_1d_2d_2( nbpac, frld , npac(1:nbpac), frld_1d (1:nbpac) , jpi, jpj ) … … 375 387 CALL tab_1d_2d_2( nbpac, rdmicif , npac(1:nbpac), rdmicif_1d(1:nbpac) , jpi, jpj ) 376 388 CALL tab_1d_2d_2( nbpac, fdvolif , npac(1:nbpac), dvlbq_1d (1:nbpac) , jpi, jpj ) 377 389 ! 378 390 ENDIF 379 391 380 392 381 ! 382 ! 393 ! Recover frld values between 0 and 1 in the Southern Hemisphere (tricky trick) 394 ! Update daily thermodynamic ice production. 383 395 !------------------------------------------------------------------------------ 384 385 396 DO jj = 1, jpj 386 397 DO ji = 1, jpi … … 392 403 IF(ln_ctl) THEN 393 404 CALL prt_ctl_info(' lim_thd end ') 394 CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif : ', tab2d_2=hsnif , clinfo2=' hsnif : ')395 CALL prt_ctl( tab2d_1=frld , clinfo1=' lim_thd: frld : ', tab2d_2=hicifp, clinfo2=' hicifp : ')396 CALL prt_ctl( tab2d_1=phicif, clinfo1=' lim_thd: phicif : ', tab2d_2=pfrld , clinfo2=' pfrld : ')397 CALL prt_ctl( tab2d_1=sist , clinfo1=' lim_thd: sist : ')398 CALL prt_ctl( tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1 : ')399 CALL prt_ctl( tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2 : ')400 CALL prt_ctl( tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3 : ')401 CALL prt_ctl( tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn : ', tab2d_2=qdtcn , clinfo2=' qdtcn : ')402 CALL prt_ctl( tab2d_1=qstoif, clinfo1=' lim_thd: qstoif : ', tab2d_2=fsbbq , clinfo2=' fsbbq : ')405 CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif : ', tab2d_2=hsnif , clinfo2=' hsnif : ' ) 406 CALL prt_ctl( tab2d_1=frld , clinfo1=' lim_thd: frld : ', tab2d_2=hicifp, clinfo2=' hicifp : ' ) 407 CALL prt_ctl( tab2d_1=phicif , clinfo1=' lim_thd: phicif : ', tab2d_2=pfrld , clinfo2=' pfrld : ' ) 408 CALL prt_ctl( tab2d_1=sist , clinfo1=' lim_thd: sist : ' ) 409 CALL prt_ctl( tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1 : ' ) 410 CALL prt_ctl( tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2 : ' ) 411 CALL prt_ctl( tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3 : ' ) 412 CALL prt_ctl( tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn : ', tab2d_2=qdtcn , clinfo2=' qdtcn : ' ) 413 CALL prt_ctl( tab2d_1=qstoif , clinfo1=' lim_thd: qstoif : ', tab2d_2=fsbbq , clinfo2=' fsbbq : ' ) 403 414 ENDIF 404 415 ! … … 422 433 & hakdif, hnzst , thth , parsub, alphs 423 434 !!------------------------------------------------------------------- 424 425 426 ! Define the initial parameters 427 ! ------------------------- 428 REWIND( numnam_ice ) 435 ! 436 REWIND( numnam_ice ) ! read namelist 429 437 READ ( numnam_ice , namicethd ) 430 IF(lwp) THEN 438 IF( lk_cpl .AND. parsub /= 0.0 ) CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' ) 439 ! 440 IF(lwp) THEN ! control print 431 441 WRITE(numout,*) 432 442 WRITE(numout,*)'lim_thd_init_2: ice parameters for ice thermodynamic computation ' … … 437 447 WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim 438 448 WRITE(numout,*)' maximum lead fraction amax = ', amax 439 WRITE(numout,*)' energy stored in brine pocket (=1) or not (=0) 449 WRITE(numout,*)' energy stored in brine pocket (=1) or not (=0) swiqst = ', swiqst 440 450 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 441 451 WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta … … 450 460 WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs 451 461 ENDIF 452 462 ! 453 463 uscomi = 1.0 / ( 1.0 - amax ) ! inverse of minimum lead fraction 454 464 rcdsn = hakdif * rcdsn 455 465 rcdic = hakdif * rcdic 456 457 IF ( ( hsndif > 100.e0 ) .OR. ( hicdif > 100.e0 )) THEN466 ! 467 IF( hsndif > 100.e0 .OR. hicdif > 100.e0 ) THEN 458 468 cnscg = 0.e0 459 469 ELSE 460 470 cnscg = rcpsn / rcpic ! ratio rcpsn/rcpic 461 471 ENDIF 462 472 ! 463 473 END SUBROUTINE lim_thd_init_2 464 474 -
trunk/NEMO/LIM_SRC_2/limthd_zdf_2.F90
r1156 r1218 22 22 USE limistate_2 23 23 USE in_out_manager 24 USE cpl_oasis3, ONLY : lk_cpl 24 25 25 26 IMPLICIT NONE … … 213 214 zghe = ( 1.0 - zihe ) * zheshth * ( 2.0 - zheshth ) & 214 215 & + zihe * 0.5 * ( 1.5 + LOG( 2.0 * zheshth ) ) 215 #if defined key_lim_cp3216 zghe = 1.0217 #endif218 216 219 217 !---effective conductivities … … 297 295 DO ji = kideb, kiut 298 296 !---computation of the derivative of energy balance function 299 #if defined key_coupled300 # if defined key_lim_cp2301 zdfts = zksndh(ji) & ! contribution of the conductive heat flux302 & + zrcpdt(ji) & ! contribution of hsu * rcp / dt303 & - dqns_ice_1d(ji) ! contribution of the total non solar radiation304 # else305 zdfts = zksndh(ji) & ! contribution of the conductive heat flux306 & + zrcpdt(ji) ! contribution of hsu * rcp / dt307 # endif308 309 #else310 297 zdfts = zksndh(ji) & ! contribution of the conductive heat flux 311 298 & + zrcpdt(ji) & ! contribution of hsu * rcp / dt 312 299 & - dqns_ice_1d (ji) ! contribution of the total non solar radiation 313 #endif314 300 !---computation of the energy balance function 315 301 zfts = - z1mi0 (ji) * qsr_ice_1d(ji) & ! net absorbed solar radiation … … 318 304 !---computation of surface temperature increment 319 305 zdts = -zfts / zdfts 320 #if defined key_lim_cp3321 zdts = zdts / 3.0322 #endif323 306 !---computation of the new surface temperature 324 307 sist_1d(ji) = sist_1d(ji) + zdts 325 326 308 END DO 327 309 … … 338 320 !---------------------------------------------------------------------- 339 321 340 DO ji = kideb, kiut 341 sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 342 #if ! defined key_coupled 343 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 344 qla_ice_1d (ji) = qla_ice_1d (ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 345 #endif 346 zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 347 END DO 322 IF ( .NOT. lk_cpl ) THEN ! duplicate the loop for performances issues 323 DO ji = kideb, kiut 324 sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 325 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 326 qla_ice_1d(ji) = qla_ice_1d(ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 327 zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 328 END DO 329 ELSE 330 DO ji = kideb, kiut 331 sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 332 zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 333 END DO 334 ENDIF 348 335 349 336 ! 5.2. Calculate available heat for surface ablation. … … 517 504 518 505 qstbif_1d(ji) = ziqf * ( qstbif_1d(ji) - zqsn_mlt_rem ) & 519 & + ( 1.0 - ziqf ) * ( qstbif_1d(ji) - qstbif_1d(ji) )506 & + ( 1.0 - ziqf ) * ( qstbif_1d(ji) - qstbif_1d(ji) ) 520 507 521 508 !-- The contribution of the energy stored in brine pockets qstbif_1d to melt … … 529 516 530 517 qstbif_1d(ji) = zihq * qstbif_1d(ji) & 531 & + ( 1.0 - zihq ) * zqstbif_old518 & + ( 1.0 - zihq ) * zqstbif_old 532 519 533 520 !--change in ice thickness due to melt at the top surface
Note: See TracChangeset
for help on using the changeset viewer.
