Changeset 3938
- Timestamp:
- 2013-06-26T09:54:16+02:00 (11 years ago)
- Location:
- branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM
- Files:
-
- 4 added
- 48 deleted
- 51 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/Job_mpi
r3935 r3938 7 7 # Type de travail 8 8 # @ job_type = parallel 9 # @ as_limit = 7.0Gb10 9 # Fichier de sortie standard 11 10 # @ output = Script_Output … … 13 12 # @ error = error 14 13 # Nombre de processus demandes 15 # @ total_tasks = 3016 # @ environment = "BATCH_NUM_PROC_TOT= 30"14 # @ total_tasks = 24 15 # @ environment = "BATCH_NUM_PROC_TOT=24" 17 16 # Temps CPU max. par processus MPI hh:mm:ss 18 17 # @ wall_clock_limit = 4:30:00 -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/xmlio_server.def
r3931 r3938 30 30 ! setting nn_nchunks_k = jpk will give a chunk size of 1 in the vertical which 31 31 ! is optimal for postprocessing which works exclusively with horizontal slabs 32 ln_nc4zip = . TRUE. ! (T) use netcdf4 chunking and compression32 ln_nc4zip = .FALSE. ! (T) use netcdf4 chunking and compression 33 33 ! (F) ignore chunking information and produce netcdf3-compatible files 34 34 / -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r2777 r3938 187 187 REAL(wp), PUBLIC :: alphaevp = 1._wp !: coeficient of the internal stresses !SB 188 188 REAL(wp), PUBLIC :: unit_fac = 1.e+09_wp !: conversion factor for ice / snow enthalpy 189 REAL(wp), PUBLIC :: hminrhg = 0.05_wp !: clem : ice thickness (in m) below which ice velocity is set to ocean velocity 189 190 190 191 ! !!** ice-salinity namelist (namicesal) ** … … 393 394 LOGICAL , PUBLIC :: ln_limdyn = .TRUE. !: flag for ice dynamics (T) or not (F) 394 395 LOGICAL , PUBLIC :: ln_nicep = .TRUE. !: flag for sea-ice points output (T) or not (F) 395 REAL(wp) , PUBLIC :: hsndif = 0.e0 !: computation of temp. in snow (0) or not (9999)396 REAL(wp) , PUBLIC :: hicdif = 0.e0 !: computation of temp. in ice (0) or not (9999)397 396 REAL(wp) , PUBLIC :: cai = 1.40e-3 !: atmospheric drag over sea ice 398 397 REAL(wp) , PUBLIC :: cao = 1.00e-3 !: atmospheric drag over ocean 399 REAL(wp) , DIMENSION(2), PUBLIC :: acrit = (/ 1.e-06 , 1.e-06 /) !: minimum fraction for leads in400 ! ! : north and south hemisphere398 REAL(wp) , PUBLIC :: amax = 0.99 !: maximum ice concentration 399 ! ! 401 400 !!-------------------------------------------------------------------------- 402 401 !! * Ice diagnostics 403 402 !!-------------------------------------------------------------------------- 404 403 !! Check if everything down here is necessary 404 LOGICAL , PUBLIC :: ln_limdiahsb = .TRUE. !: flag for ice diag (T) or not (F) 405 405 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_newice !: volume of ice formed in the leads 406 406 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dv_dt_thd !: thermodynamic growth rates … … 412 412 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_bot_me ! vertical bottom melt 413 413 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_sur_me ! vertical surface melt 414 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_res_pr ! production (growth+melt) due to limupdate 415 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vi ! transport of ice volume 414 416 INTEGER , PUBLIC :: jiindx, jjindx !: indexes of the debugging point 415 417 … … 460 462 461 463 ii = ii + 1 462 ALLOCATE( patho_case(jpi, jpj, jpl), STAT=ierr(ii) )464 ALLOCATE( patho_case(jpi,jpj,jpl) , STAT=ierr(ii) ) 463 465 464 466 ! * Ice global state variables … … 524 526 & izero (jpi,jpj,jpl) , diag_bot_gr(jpi,jpj) , diag_dyn_gr(jpi,jpj) , & 525 527 & fstroc (jpi,jpj,jpl) , diag_bot_me(jpi,jpj) , diag_sur_me(jpi,jpj) , & 526 & fhbricat (jpi,jpj,jpl) , v_newice (jpi,jpj), STAT=ierr(ii) )528 & fhbricat (jpi,jpj,jpl) , diag_res_pr(jpi,jpj) , diag_trp_vi(jpi,jpj) , v_newice(jpi,jpj) , STAT=ierr(ii) ) 527 529 528 530 ice_alloc = MAXVAL( ierr(:) ) -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
r3294 r3938 124 124 !! ** input : Namelist namicerun 125 125 !!------------------------------------------------------------------- 126 NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, a crit, hsndif, hicdif, cai, cao, ln_nicep126 NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb 127 127 !!------------------------------------------------------------------- 128 128 ! … … 140 140 WRITE(numout,*) ' ~~~~~~' 141 141 WRITE(numout,*) ' switch for ice dynamics (1) or not (0) ln_limdyn = ', ln_limdyn 142 WRITE(numout,*) ' minimum fraction for leads in the NH (SH) acrit(1/2) = ', acrit(:) 143 WRITE(numout,*) ' computation of temp. in snow (=0) or not (=9999) hsndif = ', hsndif 144 WRITE(numout,*) ' computation of temp. in ice (=0) or not (=9999) hicdif = ', hicdif 142 WRITE(numout,*) ' maximum ice concentration = ', amax 145 143 WRITE(numout,*) ' atmospheric drag over sea ice = ', cai 146 144 WRITE(numout,*) ' atmospheric drag over ocean = ', cao 147 145 WRITE(numout,*) ' Several ice points in the ice or not in ocean.output = ', ln_nicep 146 WRITE(numout,*) ' Diagnose heat/salt budget or not ln_limdiahsb = ', ln_limdiahsb 148 147 ENDIF 149 148 ! -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
r3294 r3938 23 23 USE lib_mpp ! MPP library 24 24 USE wrk_nemo ! work arrays 25 USE lib_fortran ! to use key_nosignedzero 25 26 26 27 IMPLICIT NONE … … 88 89 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 89 90 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) ) ) 90 zin0 = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj) ! Case of empty boxes & Apply mask91 zin0 = ( 1.0 - MAX( rzero, SIGN ( rone, -zslpmax ) ) ) * tms(ji,jj) ! Case of empty boxes & Apply mask 91 92 92 93 ps0 (ji,jj) = zslpmax … … 273 274 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 274 275 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) ) ) 275 zin0 = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj) ! Case of empty boxes & Apply mask276 zin0 = ( 1.0 - MAX( rzero, SIGN ( rone, -zslpmax ) ) ) * tms(ji,jj) ! Case of empty boxes & Apply mask 276 277 ! 277 278 ps0 (ji,jj) = zslpmax -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90
r2715 r3938 23 23 USE in_out_manager ! I/O manager 24 24 USE lib_mpp ! MPP library 25 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 25 26 26 27 IMPLICIT NONE … … 43 44 INTEGER :: naveg ! number of step for accumulation before averaging 44 45 REAL(wp) :: epsi06 = 1.e-6_wp ! small number 46 REAL(wp) :: epsi20 = 1.e-20_wp ! small number 45 47 46 48 CHARACTER(len= 8) :: fmtinf = '1PE13.5 ' ! format of the output values … … 70 72 !! the temporal evolution of some key variables 71 73 !!------------------------------------------------------------------- 72 INTEGER :: jv, ji, jj, jl ! dummy loop indices 73 REAL(wp) :: zshift_date ! date from the minimum ice extent 74 REAL(wp) :: zday, zday_min ! current day, day of minimum extent 75 REAL(wp) :: zafy, zamy ! temporary area of fy and my ice 74 INTEGER :: jv, ji, jj, jl ! dummy loop indices 75 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integer 76 REAL(wp) :: zshift_date ! date from the minimum ice extent 77 REAL(wp) :: zday, zday_min ! current day, day of minimum extent 78 REAL(wp) :: zafy, zamy ! temporary area of fy and my ice 76 79 REAL(wp) :: zindb 77 REAL(wp), DIMENSION(jpinfmx) :: vinfor ! temporary workingspace80 REAL(wp), DIMENSION(jpinfmx) :: vinfor ! 1D workspace 78 81 !!------------------------------------------------------------------- 79 82 80 83 ! 0) date from the minimum of ice extent 81 84 !--------------------------------------- 85 !RETURN ! use this for debugging 82 86 zday_min = 273._wp ! zday_min = date of minimum extent, here September 30th 83 87 zday = REAL(numit-nit000,wp) * rdt_ice / ( 86400._wp * REAL(nn_fsbc,wp) ) … … 112 116 ! the computation of this diagnostic is not reliable 113 117 vinfor(31) = vinfor(31) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 114 v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj) /1.0e12118 v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj) * 1.e-12 115 119 vinfor(53) = vinfor(53) + emps(ji,jj)*aire(ji,jj) * 1.e-12_wp !salt flux 116 120 vinfor(55) = vinfor(55) + fsbri(ji,jj)*aire(ji,jj) * 1.e-12_wp !brine drainage flux … … 153 157 vinfor(79) = vinfor(79) / MAX(vinfor(5),epsi06) ! 154 158 155 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(3) )) !159 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(3)+epsi20)) ! 156 160 vinfor(59) = zindb*vinfor(59) / MAX(vinfor(3),epsi06) ! divide by ice area 157 161 vinfor(61) = zindb*vinfor(61) / MAX(vinfor(3),epsi06) ! 158 162 159 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(9) )) !163 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(9)+epsi20)) ! 160 164 vinfor(65) = zindb*vinfor(65) / MAX(vinfor(9),epsi06) ! divide it by snow volume 161 165 … … 226 230 END DO 227 231 END DO 228 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(25) )) !=0 if no multiyear ice 1 if yes232 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(25)+epsi20)) !=0 if no multiyear ice 1 if yes 229 233 vinfor(49) = zindb*vinfor(49) / MAX(vinfor(25),epsi06) 230 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(27) )) !=0 if no multiyear ice 1 if yes234 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(27)+epsi20)) !=0 if no multiyear ice 1 if yes 231 235 vinfor(51) = zindb*vinfor(51) / MAX(vinfor(27),epsi06) 232 236 233 !! Fram Strait Export 234 !! 83 = area export 235 !! 84 = volume export 236 !! Fram strait in ORCA2 = 5 points 237 !! export = -v_ice*e1t*ddtb*at_i or -v_ice*e1t*ddtb*at_i*h_i 238 jj = 136 ! C grid 239 vinfor(83) = 0.0 240 vinfor(84) = 0.0 241 DO ji = 134, 138 242 vinfor(83) = vinfor(83) - v_ice(ji,jj) * & 243 e1t(ji,jj)*at_i(ji,jj)*rdt_ice * 1.e-12_wp 244 vinfor(84) = vinfor(84) - v_ice(ji,jj) * & 245 e1t(ji,jj)*vt_i(ji,jj)*rdt_ice * 1.e-12_wp 246 END DO 237 IF( cp_cfg == "orca" ) THEN !* ORCA configuration : Fram Strait Export 238 SELECT CASE ( jp_cfg ) 239 CASE ( 2 ) ! ORCA_R2 240 ij0 = 136 ; ij1 = 136 ! Fram strait : 83 = area export 241 ii0 = 134 ; ii1 = 138 ! 84 = volume export 242 DO jj = mj0(ij0),mj1(ij1) 243 DO ji = mi0(ii0),mi1(ii1) 244 vinfor(83) = vinfor(83) - v_ice(ji,jj) * e1t(ji,jj)*at_i(ji,jj)*rdt_ice * 1.e-12_wp 245 vinfor(84) = vinfor(84) - v_ice(ji,jj) * e1t(ji,jj)*vt_i(ji,jj)*rdt_ice * 1.e-12_wp 246 END DO 247 END DO 248 END SELECT 249 !!gm just above, this is NOT the correct way of evaluating the transport ! 250 !!gm mass of snow is missing and v_ice should be the mean between jj and jj+1 251 !!gm Other ORCA configurations should be added 252 ENDIF 247 253 248 254 !!------------------------------------------------------------------- … … 264 270 vinfor(32) = vinfor(32) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 265 271 v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel 272 !!gm error?? multiplication by at_i seem wrong here.... 266 273 vinfor(54) = vinfor(54) + at_i(ji,jj)*emps(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Total salt flux 267 274 vinfor(56) = vinfor(56) + at_i(ji,jj)*fsbri(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Brine drainage salt flux 268 275 vinfor(58) = vinfor(58) + at_i(ji,jj)*fseqv(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Equivalent salt flux 276 !!gm end 269 277 vinfor(60) = vinfor(60) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !SST 270 278 vinfor(62) = vinfor(62) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !SSS … … 292 300 vinfor(14) = 0.0 293 301 294 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(8) ))302 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(8)+epsi20)) 295 303 vinfor(16) = zindb * vinfor(16) / MAX(vinfor(8),epsi06) ! these have to be divided by ice vol 296 304 vinfor(30) = zindb * vinfor(30) / MAX(vinfor(8),epsi06) ! … … 298 306 vinfor(68) = zindb * vinfor(68) / MAX(vinfor(8),epsi06) ! 299 307 300 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(6) ))308 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(6)+epsi20)) 301 309 vinfor(54) = zindb * vinfor(54) / MAX(vinfor(6),epsi06) ! these have to be divided by ice extt 302 310 vinfor(56) = zindb * vinfor(56) / MAX(vinfor(6),epsi06) ! … … 305 313 ! vinfor(84) = vinfor(84) / vinfor(6) ! 306 314 307 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4) )) !315 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4)+epsi20)) ! 308 316 vinfor(60) = zindb*vinfor(60) / ( MAX(vinfor(4), epsi06) ) ! divide by ice area 309 317 vinfor(62) = zindb*vinfor(62) / ( MAX(vinfor(4), epsi06) ) ! 310 318 311 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(10) )) !319 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(10)+epsi20)) ! 312 320 vinfor(66) = zindb*vinfor(66) / MAX(vinfor(10),epsi06) ! divide it by snow volume 313 321 … … 345 353 END DO 346 354 END DO 347 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4) )) !355 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4)+epsi20)) ! 348 356 vinfor(64) = zindb * vinfor(64) / MAX(vinfor(4),epsi06) ! divide by ice extt 349 357 !! 2.2) Diagnostics dependent on age … … 377 385 END DO ! jj 378 386 END DO ! ji 379 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(26) )) !=0 if no multiyear ice 1 if yes387 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(26)+epsi20)) !=0 if no multiyear ice 1 if yes 380 388 vinfor(50) = zindb*vinfor(50) / MAX(vinfor(26),epsi06) 381 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(28) )) !=0 if no multiyear ice 1 if yes389 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(28)+epsi20)) !=0 if no multiyear ice 1 if yes 382 390 vinfor(52) = zindb*vinfor(52) / MAX(vinfor(28),epsi06) 383 391 -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r3294 r3938 28 28 USE in_out_manager ! I/O manager 29 29 USE prtctl ! Print control 30 USE lib_fortran ! glob_sum 30 31 31 32 IMPLICIT NONE … … 64 65 REAL(wp), POINTER, DIMENSION(:) :: zmsk ! i-averaged of tmask 65 66 REAL(wp), POINTER, DIMENSION(:,:) :: zu_io, zv_io ! ice-ocean velocity 66 !!--------------------------------------------------------------------- 67 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 68 !!--------------------------------------------------------------------- 67 69 68 70 CALL wrk_alloc( jpi, jpj, zu_io, zv_io ) 69 71 CALL wrk_alloc( jpj, zind, zmsk ) 72 73 ! ------------------------------- 74 !- check conservation (C Rousset) 75 IF (ln_limdiahsb) THEN 76 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 77 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 78 zchk_fw_b = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) 79 zchk_fs_b = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) 80 ENDIF 81 !- check conservation (C Rousset) 82 ! ------------------------------- 70 83 71 84 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) … … 207 220 ENDIF 208 221 ! 222 223 ! ------------------------------- 224 !- check conservation (C Rousset) 225 IF (ln_limdiahsb) THEN 226 !INTEGER :: numhsb 227 !CHARACTER (len=32) :: cl_name ! output file name 228 !cl_name = 'heat_salt_volume_budgets.txt' ! name of output file 229 230 zchk_fs = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 231 zchk_fw = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 232 233 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 234 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 235 236 IF(lwp) THEN 237 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limdyn) = ',(zchk_v_i * 86400.) 238 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * 86400.) 239 IF ( MINVAL( v_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limdyn) = ',(MINVAL(v_i) * 1.e-3) 240 IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+1.e-10 ) WRITE(numout,*) 'violation a_i>amax (limdyn) = ',MAXVAL(SUM(a_i,dim=3)) 241 ENDIF 242 !CALL ctl_opn( numhsb , cl_name , 'UNKNOWN' , 'FORMATTED' , 'SEQUENTIAL' , 1 , numout , lwp , 1 ) 243 ! 244 !WRITE( numhsb, 9010 ) "kt | heat content budget | salt content budget ", & 245 ! & "| volume budget (ssh) ", & 246 ! & "| volume budget (e3t) " 247 !WRITE( numhsb, 9010 ) " | [C] [W/m2] | [psu] [mmm/s] [SV] ", & 248 ! & "| [m3] [mmm/s] [SV] ", & 249 ! & "| [m3] [mmm/s] [SV] " 250 !IF ( kt == nitend ) CLOSE( numhsb ) 251 252 !9010 FORMAT(A80,A45,A45) 253 ENDIF 254 !- check conservation (C Rousset) 255 ! ------------------------------- 256 209 257 CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 210 258 CALL wrk_dealloc( jpj, zind, zmsk ) … … 228 276 & dm, nbiter, nbitdr, om, resl, cw, angvg, pstar, & 229 277 & c_rhg, etamn, creepl, ecc, ahi0, & 230 & nevp, telast, alphaevp 278 & nevp, telast, alphaevp, hminrhg 231 279 !!------------------------------------------------------------------- 232 280 … … 256 304 WRITE(numout,*) ' timescale for elastic waves telast = ', telast 257 305 WRITE(numout,*) ' coefficient for the solution of int. stresses alphaevp = ', alphaevp 306 WRITE(numout,*) ' min ice thickness for rheology calculations hminrhg = ', hminrhg 258 307 ENDIF 259 308 ! -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r3294 r3938 136 136 END DO ! end of sub-time step loop 137 137 138 ! ----------------------- 139 !!! final step (clem) !!! 140 DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction 141 DO ji = 1 , fs_jpim1 ! vector opt. 142 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 143 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 144 END DO 145 END DO 146 ! 147 DO jj= 2, jpjm1 ! diffusive trend : divergence of the fluxes 148 DO ji = fs_2 , fs_jpim1 ! vector opt. 149 zdiv (ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj ) & 150 & + zflv(ji,jj) - zflv(ji ,jj-1) ) / ( e1t (ji,jj) * e2t (ji,jj) ) 151 ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) ) 152 END DO 153 END DO 154 CALL lbc_lnk( ptab, 'T', 1. ) ! lateral boundary condition 155 !!! final step (clem) !!! 156 ! ----------------------- 157 138 158 IF(ln_ctl) THEN 139 159 zrlx(:,:) = ptab(:,:) - ztab0(:,:) -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r3349 r3938 5 5 !!====================================================================== 6 6 !! History : 2.0 ! 2004-01 (C. Ethe, G. Madec) Original code 7 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 7 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 8 !! - ! 2012 (C. Rousset) clean + add par_oce (for jp_sal)...bug? 8 9 !!---------------------------------------------------------------------- 9 10 #if defined key_lim3 … … 21 22 USE ice ! sea-ice variables 22 23 USE par_ice ! ice parameters 24 USE par_oce ! ocean parameters 23 25 USE dom_ice ! sea-ice domain 24 26 USE in_out_manager ! I/O manager … … 63 65 !! or from arbitrary sea-ice conditions 64 66 !!------------------------------------------------------------------- 65 INTEGER :: ji, jj, jk, jl ! dummy loop indices 66 REAL(wp) :: zeps6, zeps, ztmelts, epsi06 ! local scalars 67 REAL(wp) :: zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs 68 REAL(wp), POINTER, DIMENSION(:) :: zgfactorn, zhin 69 REAL(wp), POINTER, DIMENSION(:) :: zgfactors, zhis 67 INTEGER :: ji, jj, jk, jl ! dummy loop indices 68 INTEGER :: i_fill, jl0, ztest_1, ztest_2, ztest_3, ztest_4, ztests 69 REAL(wp) :: zarg, zV, zconv 70 REAL(wp) :: zeps06, zeps, ztmelts ! local scalars 71 REAL(wp), DIMENSION(jpl) :: zain, zhtin, zhtsn, zais, zhtis, zhtss 72 REAL(wp), DIMENSION(jpl) :: zai, zhti, zhts 70 73 REAL(wp), POINTER, DIMENSION(:,:) :: zidto ! ice indicator 71 74 !-------------------------------------------------------------------- 72 75 73 CALL wrk_alloc( jpm, zgfactorn, zgfactors, zhin, zhis )74 76 CALL wrk_alloc( jpi, jpj, zidto ) 75 77 … … 77 79 ! 1) Preliminary things 78 80 !-------------------------------------------------------------------- 79 epsi06 = 1.e-6_wp80 81 81 82 CALL lim_istate_init ! reading the initials parameters of the ice … … 106 107 107 108 ! constants for heat contents 108 zeps = 1.e-20_wp 109 zeps6 = 1.e-06_wp 110 111 ! zgfactor for initial ice distribution 112 zgfactorn(:) = 0._wp 113 zgfactors(:) = 0._wp 114 115 ! first ice type 116 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 117 zhin (1) = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 118 zgfactorn(1) = zgfactorn(1) + exp(-(zhin(1)-hginn_u)*(zhin(1)-hginn_u) * 0.5_wp ) 119 zhis (1) = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 120 zgfactors(1) = zgfactors(1) + exp(-(zhis(1)-hgins_u)*(zhis(1)-hgins_u) * 0.5_wp ) 121 END DO ! jl 122 zgfactorn(1) = aginn_u / zgfactorn(1) 123 zgfactors(1) = agins_u / zgfactors(1) 124 125 ! ------------- 126 ! new distribution, polynom of second order, conserving area and volume 127 zh1 = 0._wp 128 zh2 = 0._wp 129 zh3 = 0._wp 109 zeps = 1.e-20_wp 110 zeps06 = 1.e-06_wp 111 112 !------------------------------------------------------------------- 113 ! 3) Distribute ice concentration and thickness into the categories 114 !------------------------------------------------------------------- 115 ! snow distribution 116 zhtsn(1:jpl) = 0.d0 117 zhtss(1:jpl) = 0.d0 130 118 DO jl = 1, jpl 131 zh = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 132 zh1 = zh1 + zh 133 zh2 = zh2 + zh * zh 134 zh3 = zh3 + zh * zh * zh 119 zhtsn(jl) = hninn 120 zhtss(jl) = hnins 135 121 END DO 136 IF(lwp) WRITE(numout,*) ' zh1 : ', zh1 137 IF(lwp) WRITE(numout,*) ' zh2 : ', zh2 138 IF(lwp) WRITE(numout,*) ' zh3 : ', zh3 139 140 zvol = aginn_u * hginn_u 141 zare = aginn_u 142 IF( jpl >= 2 ) THEN 143 zbn = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3) 144 zan = ( zare - zbn*zh1 ) / zh2 122 123 ! ------------------- 124 ! Northern Hemisphere 125 ! ------------------- 126 ! initialisation of tests 127 ztest_1 = 0 128 ztest_2 = 0 129 ztest_3 = 0 130 ztest_4 = 0 131 ztests = 0 132 133 i_fill = jpl + 1 !==================================== 134 DO WHILE ( ( ztests /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 135 ! iteration !==================================== 136 i_fill = i_fill - 1 137 138 ! initialisation of ice variables for each try 139 zhtin(1:jpl) = 0.d0 140 zain (1:jpl) = 0.d0 141 142 ! *** case very thin ice: fill only category 1 143 IF ( i_fill == 1 ) THEN 144 zhtin(1) = hginn_u 145 zain (1) = aginn_u 146 147 ! *** case ice is thicker: fill categories >1 148 ELSE 149 150 ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2 151 DO jl = 1, i_fill - 1 152 zhtin(jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2. 153 END DO 154 155 ! find which category (jl0) the input ice thickness falls into 156 jl0 = i_fill 157 DO jl = 1, i_fill 158 IF ( ( hginn_u >= hi_max(jl-1) ) .AND. ( hginn_u < hi_max(jl) ) ) THEN 159 jl0 = jl 160 CYCLE 161 ENDIF 162 END DO 163 164 ! Concentrations in the (i_fill-1) categories 165 zain(jl0) = aginn_u / SQRT(REAL(jpl)) 166 DO jl = 1, i_fill - 1 167 IF ( jl == jl0 ) CYCLE 168 zarg = ( zhtin(jl) - hginn_u ) / ( hginn_u / 2. ) 169 zain(jl) = zain (jl0) * EXP(-zarg**2) 170 END DO 171 172 ! Concentration in the last (i_fill) category 173 zain(i_fill) = aginn_u - SUM( zain(1:i_fill-1) ) 174 175 ! Ice thickness in the last (i_fill) category 176 zV = SUM( zain(1:i_fill-1) * zhtin(1:i_fill-1) ) 177 zhtin(i_fill) = ( hginn_u*aginn_u - zV ) / zain(i_fill) 178 179 ENDIF ! case ice is thick or thin 180 181 !--------------------- 182 ! Compatibility tests 183 !--------------------- 184 ! Test 1: area conservation 185 zconv = ABS( aginn_u - SUM( zain(1:jpl) ) ) 186 IF ( zconv < zeps06 ) ztest_1 = 1 187 188 ! Test 2: volume conservation 189 zconv = ABS( hginn_u*aginn_u - SUM( zain(1:jpl)*zhtin(1:jpl) ) ) 190 IF ( zconv < zeps06 ) ztest_2 = 1 191 192 ! Test 3: thickness of the last category is in-bounds ? 193 IF ( zhtin(i_fill) >= hi_max(i_fill-1) ) ztest_3 = 1 194 195 ! Test 4: positivity of ice concentrations 196 ztest_4 = 1 197 DO jl = 1, i_fill 198 IF ( zain(jl) < 0.0d0 ) ztest_4 = 0 199 END DO 200 201 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 202 !============================ 203 END DO ! end iteration on categories 204 !============================ 205 ! Check if tests have passed (i.e. volume conservation...) 206 IF ( ztests .NE. 4 ) THEN 207 WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 208 WRITE(numout,*) ' !! ALERT categories distribution !!' 209 WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 210 WRITE(numout,*) ' *** ztests is not equal to 4 ' 211 WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4 145 212 ENDIF 146 213 147 IF(lwp) WRITE(numout,*) ' zvol: ', zvol 148 IF(lwp) WRITE(numout,*) ' zare: ', zare 149 IF(lwp) WRITE(numout,*) ' zbn : ', zbn 150 IF(lwp) WRITE(numout,*) ' zan : ', zan 151 152 zvol = agins_u * hgins_u 153 zare = agins_u 154 IF( jpl >= 2 ) THEN 155 zbs = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3) 156 zas = ( zare - zbs*zh1 ) / zh2 214 ! ------------------- 215 ! Southern Hemisphere 216 ! ------------------- 217 ! initialisation of tests 218 ztest_1 = 0 219 ztest_2 = 0 220 ztest_3 = 0 221 ztest_4 = 0 222 ztests = 0 223 224 i_fill = jpl + 1 !==================================== 225 DO WHILE ( ( ztests /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 226 ! iteration !==================================== 227 i_fill = i_fill - 1 228 229 ! initialisation of ice variables for each try 230 zhtis(1:jpl) = 0.d0 231 zais (1:jpl) = 0.d0 232 233 ! *** case very thin ice: fill only category 1 234 IF ( i_fill == 1 ) THEN 235 zhtis(1) = hgins_u 236 zais (1) = agins_u 237 238 ! *** case ice is thicker: fill categories >1 239 ELSE 240 241 ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2 242 DO jl = 1, i_fill - 1 243 zhtis(jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2. 244 END DO 245 246 ! find which category (jl0) the input ice thickness falls into 247 jl0 = i_fill 248 DO jl = 1, i_fill 249 IF ( ( hgins_u >= hi_max(jl-1) ) .AND. ( hgins_u < hi_max(jl) ) ) THEN 250 jl0 = jl 251 CYCLE 252 ENDIF 253 END DO 254 255 ! Concentrations in the (i_fill-1) categories 256 zais(jl0) = agins_u / SQRT(REAL(jpl)) 257 DO jl = 1, i_fill - 1 258 IF ( jl == jl0 ) CYCLE 259 zarg = ( zhtis(jl) - hgins_u ) / ( hgins_u / 2. ) 260 zais(jl) = zais (jl0) * EXP(-zarg**2) 261 END DO 262 263 ! Concentration in the last (i_fill) category 264 zais(i_fill) = agins_u - SUM( zais(1:i_fill-1) ) 265 266 ! Ice thickness in the last (i_fill) category 267 zV = SUM( zais(1:i_fill-1) * zhtis(1:i_fill-1) ) 268 zhtis(i_fill) = ( hgins_u*agins_u - zV ) / zais(i_fill) 269 270 ENDIF ! case ice is thick or thin 271 272 !--------------------- 273 ! Compatibility tests 274 !--------------------- 275 ! Test 1: area conservation 276 zconv = ABS( agins_u - SUM( zais(1:jpl) ) ) 277 IF ( zconv < zeps06 ) ztest_1 = 1 278 279 ! Test 2: volume conservation 280 zconv = ABS( hgins_u*agins_u - SUM( zais(1:jpl)*zhtis(1:jpl) ) ) 281 IF ( zconv < zeps06 ) ztest_2 = 1 282 283 ! Test 3: thickness of the last category is in-bounds ? 284 IF ( zhtis(i_fill) >= hi_max(i_fill-1) ) ztest_3 = 1 285 286 ! Test 4: positivity of ice concentrations 287 ztest_4 = 1 288 DO jl = 1, i_fill 289 IF ( zais(jl) < 0.0d0 ) ztest_4 = 0 290 END DO 291 292 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 293 !============================ 294 END DO ! end iteration on categories 295 !============================ 296 ! Check if tests have passed (i.e. volume conservation...) 297 IF ( ztests .NE. 4 ) THEN 298 WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 299 WRITE(numout,*) ' !! ALERT categories distribution !!' 300 WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 301 WRITE(numout,*) ' *** ztests is not equal to 4 ' 302 WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4 157 303 ENDIF 158 159 IF(lwp) WRITE(numout,*) ' zvol: ', zvol 160 IF(lwp) WRITE(numout,*) ' zare: ', zare 161 IF(lwp) WRITE(numout,*) ' zbn : ', zbn 162 IF(lwp) WRITE(numout,*) ' zan : ', zan 163 164 !end of new lines 165 ! ------------- 166 !!! 167 ! retour a LIMA_MEC 168 ! ! second ice type 169 ! zdummy = hi_max(ice_cat_bounds(2,1)-1) 170 ! hi_max(ice_cat_bounds(2,1)-1) = 0.0 171 172 ! ! here to change !!!! 173 ! jm = 2 174 ! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 175 ! zhin (2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 176 ! zhin (2) = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm ) + & 177 ! hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm) ) / 2.0 178 ! zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0) 179 ! zhis (2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 180 ! zhis (2) = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm ) + & 181 ! hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm) ) / 2.0 182 ! zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0) 183 ! END DO ! jl 184 ! zgfactorn(2) = aginn_d / zgfactorn(2) 185 ! zgfactors(2) = agins_d / zgfactors(2) 186 ! hi_max(ice_cat_bounds(2,1)-1) = zdummy 187 ! END retour a LIMA_MEC 188 !!! 189 190 !!gm optimisation : loop over the ice categories inside the ji, jj loop !!! 191 304 !! 305 !! 306 !! 307 ! --------------------- 308 ! 4) fill ice variables 309 ! --------------------- 310 zai(:) = 0.0 311 zhti(:)= 0.0 312 zhts(:)= 0.0 192 313 DO jj = 1, jpj 193 314 DO ji = 1, jpi … … 195 316 !--- Northern hemisphere 196 317 !---------------------------------------------------------------- 197 IF( fcor(ji,jj) >= 0._wp ) THEN 198 318 IF( fcor(ji,jj) >= 0._wp ) THEN 319 zai(:) = zain(:) ; zhti(:) = zhtin(:) ; zhts(:) = zhtsn(:) 320 ELSE 321 zai(:) = zais(:) ; zhti(:) = zhtis(:) ; zhts(:) = zhtss(:) 322 ENDIF 323 324 DO jl = 1, jpl 199 325 !----------------------- 200 326 ! Ice area / thickness 201 327 !----------------------- 202 203 IF ( jpl .EQ. 1) THEN ! one category 204 205 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 206 a_i(ji,jj,jl) = zidto(ji,jj) * aginn_u 207 ht_i(ji,jj,jl) = zidto(ji,jj) * hginn_u 208 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 209 END DO 210 211 ELSE ! several categories 212 213 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 214 zhin(1) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 215 a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* & 216 (zhin(1)-hginn_u)/2.0) , epsi06) 217 ! new line 218 a_i(ji,jj,jl) = zidto(ji,jj) * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) ) 219 ht_i(ji,jj,jl) = zidto(ji,jj) * zhin(1) 220 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 221 END DO 222 223 ENDIF 224 225 226 !!! 227 ! retour a LIMA_MEC 228 ! !ridged ice 229 ! zdummy = hi_max(ice_cat_bounds(2,1)-1) 230 ! hi_max(ice_cat_bounds(2,1)-1) = 0.0 231 ! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories 232 ! zhin(2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 233 ! a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* & 234 ! (zhin(2)-hginn_d)/2.0) , epsi06) 235 ! ht_i(ji,jj,jl) = zidto(ji,jj) * zhin(2) 236 ! v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 237 ! END DO 238 ! hi_max(ice_cat_bounds(2,1)-1) = zdummy 239 240 ! !rafted ice 241 ! jl = 6 242 ! a_i(ji,jj,jl) = 0.0 243 ! ht_i(ji,jj,jl) = 0.0 244 ! v_i(ji,jj,jl) = 0.0 245 ! END retour a LIMA_MEC 246 !!! 247 248 DO jl = 1, jpl 249 250 !------------- 251 ! Snow depth 252 !------------- 253 ht_s(ji,jj,jl) = zidto(ji,jj) * hninn 254 v_s(ji,jj,jl) = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 255 256 !--------------- 257 ! Ice salinity 258 !--------------- 259 sm_i(ji,jj,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1 260 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 261 262 !---------- 263 ! Ice age 264 !---------- 265 o_i(ji,jj,jl) = zidto(ji,jj) * 1.0 + ( 1.0 - zidto(ji,jj) ) 266 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) 267 268 !------------------------------ 269 ! Sea ice surface temperature 270 !------------------------------ 271 272 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 273 274 !------------------------------------ 275 ! Snow temperature and heat content 276 !------------------------------------ 277 278 DO jk = 1, nlay_s 279 t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 280 ! Snow energy of melting 281 e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 282 ! Change dimensions 283 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 284 ! Multiply by volume, so that heat content in 10^9 Joules 285 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 286 v_s(ji,jj,jl) / nlay_s 287 END DO !jk 288 289 !----------------------------------------------- 290 ! Ice salinities, temperature and heat content 291 !----------------------------------------------- 292 293 DO jk = 1, nlay_i 294 t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 295 s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1 296 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 297 298 ! heat content per unit volume 299 e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 300 ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 301 + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 302 - rcp * ( ztmelts - rtt ) & 303 ) 304 305 ! Correct dimensions to avoid big values 306 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 307 308 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 309 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 310 area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 311 nlay_i 312 END DO ! jk 313 314 END DO ! jl 315 316 ELSE ! on fcor 317 318 !--- Southern hemisphere 319 !---------------------------------------------------------------- 320 321 !----------------------- 322 ! Ice area / thickness 323 !----------------------- 324 325 IF ( jpl .EQ. 1) THEN ! one category 326 327 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 328 a_i(ji,jj,jl) = zidto(ji,jj) * agins_u 329 ht_i(ji,jj,jl) = zidto(ji,jj) * hgins_u 330 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 331 END DO 332 333 ELSE ! several categories 334 335 !level ice 336 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories 337 338 zhis(1) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 339 a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * & 340 (zhis(1)-hgins_u)/2.0) , epsi06 ) 341 ! new line square distribution volume conserving 342 a_i(ji,jj,jl) = zidto(ji,jj) * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) ) 343 ht_i(ji,jj,jl) = zidto(ji,jj) * zhis(1) 344 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 345 346 END DO ! jl 347 348 ENDIF 349 350 !!! 351 ! retour a LIMA_MEC 352 ! !ridged ice 353 ! zdummy = hi_max(ice_cat_bounds(2,1)-1) 354 ! hi_max(ice_cat_bounds(2,1)-1) = 0.0 355 ! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories 356 ! zhis(2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 357 ! a_i(ji,jj,jl) = zidto(ji,jj)*MAX( zgfactors(2) & 358 ! & * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 ) 359 ! ht_i(ji,jj,jl) = zidto(ji,jj) * zhis(2) 360 ! v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 361 ! END DO 362 ! hi_max(ice_cat_bounds(2,1)-1) = zdummy 363 364 ! !rafted ice 365 ! jl = 6 366 ! a_i(ji,jj,jl) = 0.0 367 ! ht_i(ji,jj,jl) = 0.0 368 ! v_i(ji,jj,jl) = 0.0 369 ! END retour a LIMA_MEC 370 !!! 371 372 DO jl = 1, jpl !over thickness categories 373 374 !--------------- 375 ! Snow depth 376 !--------------- 377 378 ht_s(ji,jj,jl) = zidto(ji,jj) * hnins 379 v_s(ji,jj,jl) = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 380 381 !--------------- 382 ! Ice salinity 383 !--------------- 384 385 sm_i(ji,jj,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1 386 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 387 388 !---------- 389 ! Ice age 390 !---------- 391 392 o_i(ji,jj,jl) = zidto(ji,jj) * 1.0 + ( 1.0 - zidto(ji,jj) ) 393 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) 394 395 !------------------------------ 396 ! Sea ice surface temperature 397 !------------------------------ 398 399 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 400 401 !---------------------------------- 402 ! Snow temperature / heat content 403 !---------------------------------- 404 405 DO jk = 1, nlay_s 406 t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 407 ! Snow energy of melting 408 e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 409 ! Change dimensions 410 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 411 ! Multiply by volume, so that heat content in 10^9 Joules 412 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 413 v_s(ji,jj,jl) / nlay_s 414 END DO 415 416 !--------------------------------------------- 417 ! Ice temperature, salinity and heat content 418 !--------------------------------------------- 419 420 DO jk = 1, nlay_i 421 t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 422 s_i(ji,jj,jk,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1 423 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 424 425 ! heat content per unit volume 426 e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 427 ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 428 + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 429 - rcp * ( ztmelts - rtt ) & 430 ) 431 432 ! Correct dimensions to avoid big values 433 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 434 435 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 436 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 437 area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 438 nlay_i 439 END DO !jk 440 441 END DO ! jl 442 443 ENDIF ! on fcor 444 328 a_i(ji,jj,jl) = zidto(ji,jj) * zai(jl) 329 ht_i(ji,jj,jl) = zidto(ji,jj) * zhti(jl) 330 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 331 332 !------------- 333 ! Snow depth 334 !------------- 335 ht_s(ji,jj,jl) = zidto(ji,jj) * zhts(jl) 336 v_s(ji,jj,jl) = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 337 338 !--------------- 339 ! Ice salinity 340 !--------------- 341 sm_i(ji,jj,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1 342 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 343 344 !---------- 345 ! Ice age 346 !---------- 347 o_i(ji,jj,jl) = zidto(ji,jj) * 1.0 + ( 1.0 - zidto(ji,jj) ) 348 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) 349 350 !------------------------------ 351 ! Sea ice surface temperature 352 !------------------------------ 353 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 354 355 !------------------------------------ 356 ! Snow temperature and heat content 357 !------------------------------------ 358 DO jk = 1, nlay_s 359 t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 360 ! Snow energy of melting 361 e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 362 ! Change dimensions 363 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 364 ! Multiply by volume, so that heat content in 10^9 Joules 365 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 366 v_s(ji,jj,jl) / REAL( nlay_s ) 367 END DO !jk 368 369 !----------------------------------------------- 370 ! Ice salinities, temperature and heat content 371 !----------------------------------------------- 372 DO jk = 1, nlay_i 373 t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 374 s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1 375 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 376 377 ! heat content per unit volume 378 e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 379 ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 380 + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 381 - rcp * ( ztmelts - rtt ) & 382 ) 383 384 ! Correct dimensions to avoid big values 385 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 386 387 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 388 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 389 area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 390 REAL( nlay_i ) 391 END DO ! jk 392 393 END DO ! jl 394 445 395 END DO 446 396 END DO 447 448 !-------------------------------------------------------------------- 449 ! 3) Global ice variables for output diagnostics | 397 398 IF(lwp) THEN ! control print 399 WRITE(numout,*) '~~~~~~~~~~~~~~~' 400 WRITE(numout,*) ' max ice thickness = ', (MAXVAL(ht_i(:,:,jl)),jl=1,jpl) 401 WRITE(numout,*) ' max snow thickness = ', (MAXVAL(ht_s(:,:,jl)),jl=1,jpl) 402 WRITE(numout,*) ' max ice area = ', (MAXVAL(a_i(:,:,jl)),jl=1,jpl) 403 WRITE(numout,*) ' max ice salinity = ', (MAXVAL(sm_i(:,:,jl)),jl=1,jpl) 404 WRITE(numout,*) ' min ice salinity = ', (MINVAL(sm_i(:,:,jl)),jl=1,jpl) 405 WRITE(numout,*) ' max ice surf temper = ', (MAXVAL(t_su(:,:,jl)-rtt),jl=1,jpl) 406 WRITE(numout,*) ' min ice surf temper = ', (MINVAL(t_su(:,:,jl)-rtt),jl=1,jpl) 407 ENDIF 408 409 !-------------------------------------------------------------------- 410 ! 5) Global ice variables for output diagnostics | 450 411 !-------------------------------------------------------------------- 451 412 … … 458 419 459 420 !-------------------------------------------------------------------- 460 ! 4) Moments for advection421 ! 6) Moments for advection 461 422 !-------------------------------------------------------------------- 462 423 … … 485 446 sxysal (:,:,:) = 0.e0 486 447 487 !-------------------------------------------------------------------- 488 ! 5) Lateral boundary conditions | 448 sxopw(:,:) = 0.e0 ; syopw(:,:) = 0.e0; sxxopw(:,:) = 0.e0; syyopw(:,:) = 0.e0; sxyopw(:,:) = 0.e0 449 450 !-------------------------------------------------------------------- 451 ! 7) Lateral boundary conditions | 489 452 !-------------------------------------------------------------------- 490 453 … … 518 481 CALL lbc_lnk( fsbbq , 'T', 1. ) 519 482 ! 520 CALL wrk_dealloc( jpm, zgfactorn, zgfactors, zhin, zhis )521 483 CALL wrk_dealloc( jpi, jpj, zidto ) 522 484 ! -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r3294 r3938 28 28 USE wrk_nemo ! work arrays 29 29 USE prtctl ! Print control 30 ! Check budget (Rousset) 31 USE iom ! I/O manager 32 USE lib_fortran ! glob_sum 33 USE limdiahsb 30 34 31 35 IMPLICIT NONE … … 50 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/ 51 55 ! ! closing associated w/ category n 52 53 56 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness 54 57 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness … … 138 141 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2) 139 142 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 143 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 144 ! mass and salt flux (clem) 145 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... 140 146 !!----------------------------------------------------------------------------- 141 147 142 148 CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 149 150 CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem 143 151 144 152 IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only) … … 148 156 CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ') 149 157 ENDIF 158 159 IF( ln_limdyn ) THEN ! Start ridging and rafting ! 160 ! ------------------------------- 161 !- check conservation (C Rousset) 162 IF (ln_limdiahsb) THEN 163 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 164 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 165 zchk_fw_b = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) 166 zchk_fs_b = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) 167 ENDIF 168 !- check conservation (C Rousset) 169 ! ------------------------------- 170 171 ! mass and salt flux init (clem) 172 zviold(:,:,:) = v_i(:,:,:) 173 zvsold(:,:,:) = v_s(:,:,:) 174 zsmvold(:,:,:) = smv_i(:,:,:) 150 175 151 176 !-----------------------------------------------------------------------------! … … 201 226 ! to give asum = 1.0 after ridging. 202 227 203 divu_adv(ji,jj) = ( 1._wp- asum(ji,jj) ) / rdt_ice ! asum found in ridgeprep228 divu_adv(ji,jj) = ( amax - asum(ji,jj) ) / rdt_ice ! asum found in ridgeprep 204 229 205 230 IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) … … 286 311 DO jj = 1, jpj 287 312 DO ji = 1, jpi 288 IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN313 IF (ABS(asum(ji,jj) - amax ) .LT. epsi11) THEN 289 314 closing_net(ji,jj) = 0._wp 290 315 opning (ji,jj) = 0._wp 291 316 ELSE 292 317 iterate_ridging = 1 293 divu_adv (ji,jj) = ( 1._wp - asum(ji,jj)) / rdt_ice318 divu_adv (ji,jj) = ( amax - asum(ji,jj) ) / rdt_ice 294 319 closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 295 320 opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) ) … … 330 355 DO ji = 1, jpi 331 356 332 IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true.357 IF (ABS(asum(ji,jj) - amax) .GT. epsi11) asum_error = .true. 333 358 334 359 dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti … … 349 374 DO jj = 1, jpj 350 375 DO ji = 1, jpi 351 IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) THEN ! there is a bug376 IF (ABS(asum(ji,jj) - amax) .GT. epsi11) THEN ! there is a bug 352 377 WRITE(numout,*) ' ' 353 378 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) … … 377 402 CALL lim_var_glo2eqv 378 403 CALL lim_itd_me_zapsmall 404 405 !-------------------------------- 406 ! Update mass/salt fluxes (clem) 407 !-------------------------------- 408 DO jl = 1, jpl 409 DO jj = 1, jpj 410 DO ji = 1, jpi 411 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) / rdt_ice 412 rdmicif(ji,jj) = rdmicif(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 413 rdmsnif(ji,jj) = rdmsnif(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 414 fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic / rdt_ice 415 END DO 416 END DO 417 END DO 379 418 380 419 !----------------- … … 425 464 ENDIF 426 465 466 ! ------------------------------- 467 !- check conservation (C Rousset) 468 IF (ln_limdiahsb) THEN 469 zchk_fs = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 470 zchk_fw = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 471 472 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 473 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 474 475 IF(lwp) THEN 476 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_me) = ',(zchk_v_i * 86400.) 477 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * 86400.) 478 IF ( MINVAL( v_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_me) = ',(MINVAL(v_i) * 1.e-3) 479 IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_me) = ',MAXVAL(SUM(a_i,dim=3)) 480 IF ( MINVAL( a_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_me) = ',MINVAL(a_i) 481 ENDIF 482 ENDIF 483 !- check conservation (C Rousset) 484 ! ------------------------------- 485 427 486 !-------------------------! 428 487 ! Back to initial values 429 488 !-------------------------! 430 431 489 ! update of fields will be made later in lim update 432 490 u_ice(:,:) = old_u_ice(:,:) … … 439 497 oa_i(:,:,:) = old_oa_i(:,:,:) 440 498 IF( num_sal == 2 .OR. num_sal == 4 ) smv_i(:,:,:) = old_smv_i(:,:,:) 441 442 499 !----------------------------------------------------! 443 500 ! Advection of ice in a free cell, newly ridged ice … … 448 505 449 506 ! heat content has to be corrected before ice volume 450 DO jl = 1, jpl 451 DO jk = 1, nlay_i 452 DO jj = 1, jpj 453 DO ji = 1, jpi 454 IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 455 ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 456 old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl) 457 d_e_i_trp(ji,jj,jk,jl) = 0._wp 458 ENDIF 459 END DO 460 END DO 461 END DO 462 END DO 463 464 DO jl = 1, jpl 465 DO jj = 1, jpj 466 DO ji = 1, jpi 467 IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 468 ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 469 old_v_i(ji,jj,jl) = d_v_i_trp(ji,jj,jl) 470 d_v_i_trp(ji,jj,jl) = 0._wp 471 old_a_i(ji,jj,jl) = d_a_i_trp(ji,jj,jl) 472 d_a_i_trp(ji,jj,jl) = 0._wp 473 old_v_s(ji,jj,jl) = d_v_s_trp(ji,jj,jl) 474 d_v_s_trp(ji,jj,jl) = 0._wp 475 old_e_s(ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 476 d_e_s_trp(ji,jj,1,jl) = 0._wp 477 old_oa_i(ji,jj,jl) = d_oa_i_trp(ji,jj,jl) 478 d_oa_i_trp(ji,jj,jl) = 0._wp 479 IF( num_sal == 2 .OR. num_sal == 4 ) old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 480 d_smv_i_trp(ji,jj,jl) = 0._wp 481 ENDIF 482 END DO 483 END DO 484 END DO 485 507 !clem@order 508 ! DO jl = 1, jpl 509 ! DO jk = 1, nlay_i 510 ! DO jj = 1, jpj 511 ! DO ji = 1, jpi 512 ! IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 513 ! ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 514 ! old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl) 515 ! d_e_i_trp(ji,jj,jk,jl) = 0._wp 516 ! ENDIF 517 ! END DO 518 ! END DO 519 ! END DO 520 ! END DO 521 ! 522 ! DO jl = 1, jpl 523 ! DO jj = 1, jpj 524 ! DO ji = 1, jpi 525 ! IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 526 ! ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 527 ! old_v_i(ji,jj,jl) = d_v_i_trp(ji,jj,jl) 528 ! d_v_i_trp(ji,jj,jl) = 0._wp 529 ! old_a_i(ji,jj,jl) = d_a_i_trp(ji,jj,jl) 530 ! d_a_i_trp(ji,jj,jl) = 0._wp 531 ! old_v_s(ji,jj,jl) = d_v_s_trp(ji,jj,jl) 532 ! d_v_s_trp(ji,jj,jl) = 0._wp 533 ! old_e_s(ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 534 ! d_e_s_trp(ji,jj,1,jl) = 0._wp 535 ! old_oa_i(ji,jj,jl) = d_oa_i_trp(ji,jj,jl) 536 ! d_oa_i_trp(ji,jj,jl) = 0._wp 537 ! IF( num_sal == 2 .OR. num_sal == 4 ) old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 538 ! d_smv_i_trp(ji,jj,jl) = 0._wp 539 ! ENDIF 540 ! END DO 541 ! END DO 542 ! END DO 543 !clem@order 544 ENDIF ! ln_limdyn=.true. 545 ! 486 546 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 547 ! 548 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem 487 549 ! 488 550 END SUBROUTINE lim_itd_me … … 1086 1148 afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 1087 1149 1088 IF (afrac(ji,jj) > 1.0+ epsi11) THEN !riging1150 IF (afrac(ji,jj) > amax + epsi11) THEN !riging 1089 1151 large_afrac = .true. 1090 ELSEIF (afrac(ji,jj) > 1.0) THEN ! roundoff error1091 afrac(ji,jj) = 1.01152 ELSEIF (afrac(ji,jj) > amax) THEN ! roundoff error 1153 afrac(ji,jj) = amax 1092 1154 ENDIF 1093 IF (afrft(ji,jj) > 1.0+ epsi11) THEN !rafting1155 IF (afrft(ji,jj) > amax + epsi11) THEN !rafting 1094 1156 large_afrft = .true. 1095 ELSEIF (afrft(ji,jj) > 1.0) THEN ! roundoff error1096 afrft(ji,jj) = 1.01157 ELSEIF (afrft(ji,jj) > amax) THEN ! roundoff error 1158 afrft(ji,jj) = amax 1097 1159 ENDIF 1098 1160 … … 1137 1199 1138 1200 ! ! excess of salt is flushed into the ocean 1139 fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic / rdt_ice 1201 !fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic / rdt_ice 1202 1203 !rdmicif(ji,jj) = rdmicif(ji,jj) + vsw(ji,jj) * rhoic ! gurvan: increase in ice volume du to seawater frozen in voids 1140 1204 1141 1205 !------------------------------------ … … 1148 1212 dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 1149 1213 dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 1150 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) / rdt_ice1214 !clem diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) / rdt_ice 1151 1215 opening (ji,jj) = opening (ji,jj) + opning(ji,jj)*rdt_ice 1152 1216 … … 1217 1281 1218 1282 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 1219 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / nlay_i1283 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 1220 1284 1221 1285 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) … … 1240 1304 ji = indxi(ij) 1241 1305 jj = indxj(ij) 1242 IF( afrac(ji,jj) > 1.0+ epsi11 ) THEN1306 IF( afrac(ji,jj) > amax + epsi11 ) THEN 1243 1307 WRITE(numout,*) '' 1244 1308 WRITE(numout,*) ' ardg > a_i' … … 1252 1316 ji = indxi(ij) 1253 1317 jj = indxj(ij) 1254 IF( afrft(ji,jj) > 1.0+ epsi11 ) THEN1318 IF( afrft(ji,jj) > amax + epsi11 ) THEN 1255 1319 WRITE(numout,*) '' 1256 1320 WRITE(numout,*) ' arft > a_i' -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r3294 r3938 34 34 USE lib_mpp ! MPP library 35 35 USE wrk_nemo ! work arrays 36 USE lib_fortran ! to use key_nosignedzero 36 37 37 38 IMPLICIT NONE … … 44 45 PUBLIC lim_itd_shiftice 45 46 46 REAL(wp) :: epsi20 = 1 e-20_wp ! constant values47 REAL(wp) :: epsi13 = 1 e-13_wp !48 REAL(wp) :: epsi10 = 1 e-10_wp !47 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values 48 REAL(wp) :: epsi13 = 1.e-13_wp ! 49 REAL(wp) :: epsi10 = 1.e-10_wp ! 49 50 50 51 !!---------------------------------------------------------------------- … … 66 67 ! 67 68 INTEGER :: jl, ja, jm, jbnd1, jbnd2 ! ice types dummy loop index 68 69 !!------------------------------------------------------------------ 69 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 70 !!------------------------------------------------------------------ 71 ! ------------------------------- 72 !- check conservation (C Rousset) 73 IF (ln_limdiahsb) THEN 74 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 75 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 76 zchk_fw_b = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) 77 zchk_fs_b = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) 78 ENDIF 79 !- check conservation (C Rousset) 80 ! ------------------------------- 70 81 71 82 IF( kt == nit000 .AND. lwp ) THEN … … 106 117 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 107 118 d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 108 119 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 109 120 d_smv_i_thd(:,:,:) = 0._wp 110 121 IF( num_sal == 2 .OR. num_sal == 4 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 122 123 ! diag only (clem) 124 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) / rdt_ice * 86400.0 111 125 112 126 IF(ln_ctl) THEN ! Control print … … 141 155 END DO 142 156 ENDIF 143 157 ! 158 ! ------------------------------- 159 !- check conservation (C Rousset) 160 IF (ln_limdiahsb) THEN 161 zchk_fs = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 162 zchk_fw = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 163 164 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 165 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 166 167 IF(lwp) THEN 168 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_th) = ',(zchk_v_i * 86400.) 169 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * 86400.) 170 IF ( MINVAL( v_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_th) = ',(MINVAL(v_i) * 1.e-3) 171 IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_th) = ',MAXVAL(SUM(a_i,dim=3)) 172 ENDIF 173 ENDIF 174 !- check conservation (C Rousset) 175 ! ------------------------------- 176 ! 144 177 !- Recover Old values 145 178 a_i(:,:,:) = old_a_i (:,:,:) … … 148 181 e_s(:,:,:,:) = old_e_s (:,:,:,:) 149 182 e_i(:,:,:,:) = old_e_i (:,:,:,:) 150 ! 183 !?? oa_i(:,:,:) = old_oa_i(:,:,:) 151 184 IF( num_sal == 2 .OR. num_sal == 4 ) smv_i(:,:,:) = old_smv_i (:,:,:) 152 185 ! … … 239 272 DO jj = 1, jpj 240 273 DO ji = 1, jpi 241 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl) )) !0 if no ice and 1 if yes274 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 242 275 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi10) * zindb 243 zindb = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl) )) !0 if no ice and 1 if yes276 zindb = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 244 277 zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),epsi10) * zindb 245 278 IF( a_i(ji,jj,jl) > 1e-6 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) … … 405 438 * a_i(zji,zjj,klbnd) / ( a_i(zji,zjj,klbnd) - zda0 ) 406 439 a_i(zji,zjj,klbnd) = a_i(zji,zjj,klbnd) - zda0 407 v_i(zji,zjj,klbnd) = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd) 440 v_i(zji,zjj,klbnd) = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd) ! clem@useless 408 441 ENDIF ! zetamax > 0 409 442 ! ji, a_i > epsi10 … … 499 532 a_i(zji,zjj,1) = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin 500 533 ht_i(zji,zjj,1) = zhimin 501 v_i(zji,zjj,1) = a_i(zji,zjj,1)*ht_i(zji,zjj,1) 534 v_i(zji,zjj,1) = a_i(zji,zjj,1)*ht_i(zji,zjj,1) !clem@useless 502 535 ENDIF 503 536 END DO !ji … … 859 892 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 860 893 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 861 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl) )) !0 if no ice and 1 if yes894 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 862 895 ELSE 863 896 ht_i(ji,jj,jl) = 0._wp … … 968 1001 zdaice(ji,jj,jl) = a_i(ji,jj,jl) 969 1002 zdvice(ji,jj,jl) = v_i(ji,jj,jl) 1003 ! begin TECLIM change 1004 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl) 1005 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl) 1006 zdaice(ji,jj,jl) = a_i(ji,jj,jl)/2 1007 zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1))/2 1008 ! end TECLIM change 970 1009 ENDIF 971 1010 END DO ! ji -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r3294 r3938 35 35 USE dom_ice_2 ! LIM2: ice domain 36 36 #endif 37 USE lib_fortran ! to use key_nosignedzero 38 39 #if defined key_bdy 40 USE bdyice_lim 41 #endif 37 42 38 43 IMPLICIT NONE … … 43 48 REAL(wp) :: rzero = 0._wp ! constant values 44 49 REAL(wp) :: rone = 1._wp ! constant values 50 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values 45 51 46 52 !! * Substitutions … … 189 195 #endif 190 196 ! tmi = 1 where there is ice or on land 191 tmi(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd) ) ) * tms(ji,jj)197 tmi(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) ) ) ) * tms(ji,jj) 192 198 END DO 193 199 END DO … … 569 575 570 576 ENDIF 577 578 !#if defined key_bdy 579 ! ! clem: change u_ice and v_ice at the boundary for each iteration 580 ! CALL bdy_ice_lim_dyn() 581 !#endif 571 582 572 583 IF(ln_ctl) THEN … … 580 591 ENDIF 581 592 582 ! 593 ! ! ==================== ! 583 594 END DO ! end loop over jter ! 584 595 ! ! ==================== ! 585 586 596 ! 587 597 !------------------------------------------------------------------------------! 588 598 ! 4) Prevent ice velocities when the ice is thin 589 599 !------------------------------------------------------------------------------! 590 ! 591 ! If the ice thickness is below 1cm then ice velocity should equal the 600 !clem : add hminrhg in the namelist 601 ! 602 ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 592 603 ! ocean velocity, 593 604 ! This prevents high velocity when ice is thin … … 598 609 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 599 610 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 600 IF ( zdummy .LE. 5.0e-2) THEN611 IF ( zdummy .LE. hminrhg ) THEN 601 612 u_ice(ji,jj) = u_oce(ji,jj) 602 613 v_ice(ji,jj) = v_oce(ji,jj) … … 607 618 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 608 619 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 620 621 ! clem: change u_ice and v_ice at the boundary 622 #if defined key_bdy 623 CALL bdy_ice_lim_dyn() 624 #endif 609 625 610 626 DO jj = k_j1+1, k_jpj-1 … … 612 628 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 613 629 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 614 IF ( zdummy .LE. 5.0e-2) THEN630 IF ( zdummy .LE. hminrhg ) THEN 615 631 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 616 632 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & … … 637 653 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 638 654 639 IF ( zdummy .LE. 5.0e-2) THEN655 IF ( zdummy .LE. hminrhg ) THEN 640 656 641 657 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & … … 692 708 divu_i (ji,jj) = zdd (ji,jj) 693 709 delta_i(ji,jj) = deltat(ji,jj) 694 shear_i(ji,jj) = zds (ji,jj) 710 ! begin TECLIM change 711 ! shear_i(ji,jj) = zds (ji,jj) 712 zdst = ( e2u( ji , jj ) * v_ice1(ji,jj) & 713 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 714 & + e1v( ji , jj ) * u_ice2(ji,jj) & 715 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) & 716 & ) & 717 & / area(ji,jj) 718 shear_i(ji,jj) = SQRT( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) 719 ! end TECLIM change 695 720 END DO 696 721 END DO … … 699 724 CALL lbc_lnk( divu_i (:,:), 'T', 1. ) 700 725 CALL lbc_lnk( delta_i(:,:), 'T', 1. ) 701 CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 726 ! CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 727 CALL lbc_lnk( shear_i(:,:), 'T', 1. ) 702 728 703 729 ! * Store the stress tensor for the next time step -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
r3294 r3938 25 25 USE lib_mpp ! MPP library 26 26 USE wrk_nemo ! work arrays 27 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 27 28 28 29 IMPLICIT NONE … … 161 162 CALL iom_rstput( iter, nitrst, numriw, 'v_ice' , v_ice ) 162 163 CALL iom_rstput( iter, nitrst, numriw, 'fsbbq' , fsbbq ) 164 CALL iom_rstput( iter, nitrst, numriw, 'iatte' , iatte ) ! clem modif 165 CALL iom_rstput( iter, nitrst, numriw, 'oatte' , oatte ) ! clem modif 163 166 CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) 164 167 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) … … 339 342 !Control of date 340 343 341 IF( ( nit000 - INT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 ) &344 IF( ( nit000 - NINT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 ) & 342 345 & CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nit000 in ice restart', & 343 346 & ' verify the file or rerun with the value 0 for the', & 344 347 & ' control of time parameter nrstdt' ) 345 IF( INT(zfice) /= nn_fsbc .AND. ABS( nrstdt ) == 1 ) &348 IF( NINT(zfice) /= nn_fsbc .AND. ABS( nrstdt ) == 1 ) & 346 349 & CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nn_fsbc in ice restart', & 347 350 & ' verify the file or rerun with the value 0 for the', & … … 436 439 CALL iom_get( numrir, jpdom_autoglo, 'v_ice' , v_ice ) 437 440 CALL iom_get( numrir, jpdom_autoglo, 'fsbbq' , fsbbq ) 441 CALL iom_get( numrir, jpdom_autoglo, 'iatte' , iatte ) ! clem modif 442 CALL iom_get( numrir, jpdom_autoglo, 'oatte' , oatte ) ! clem modif 438 443 CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i ) 439 444 CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i ) … … 562 567 END DO 563 568 ! 564 CALL iom_close( numrir )569 !clem CALL iom_close( numrir ) 565 570 ! 566 571 CALL wrk_dealloc( nlay_i, zs_zero ) -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r3294 r3938 10 10 !! ! + simplification of the ice-ocean stress calculation 11 11 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 12 !! - ! 2012 (D. Iovino) salt flux change 13 !! - ! 2012-05 (C. Rousset) add penetration solar flux 12 14 !!---------------------------------------------------------------------- 13 15 #if defined key_lim3 … … 34 36 USE prtctl ! Print control 35 37 USE cpl_oasis3, ONLY : lk_cpl 38 USE traqsr ! clem: add penetration of solar flux into the calculation of heat budget 39 USE lib_fortran ! to use key_nosignedzero 36 40 37 41 IMPLICIT NONE … … 44 48 REAL(wp) :: r1_rdtice ! = 1. / rdt_ice 45 49 REAL(wp) :: epsi16 = 1.e-16_wp ! constant values 50 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values 46 51 REAL(wp) :: rzero = 0._wp 47 52 REAL(wp) :: rone = 1._wp … … 100 105 INTEGER :: ifvt, i1mfr, idfr ! some switches 101 106 INTEGER :: iflt, ial, iadv, ifral, ifrdv 102 REAL(wp) :: zinda, zfons, zpme ! local scalars 107 REAL(wp) :: zinda, zindb, zfons, zpme ! local scalars 108 REAL(wp) :: zfmm ! IOVINO freezing minus melting (F-M) 103 109 REAL(wp), POINTER, DIMENSION(:,:) :: zfcm1 , zfcm2 ! solar/non solar heat fluxes 104 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb, zalbp ! 2D/3D workspace … … 117 123 DO ji = 1, jpi 118 124 zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 119 ifvt = zinda * MAX( rzero , SIGN( rone, -phicif (ji,jj) ) ) !subscripts are bad here 120 i1mfr = 1.0 - MAX( rzero , SIGN( rone , - ( at_i(ji,jj) ) ) ) 125 zindb = 1.0 - MAX( rzero , SIGN( rone , - iatte(ji,jj) ) ) 126 ifvt = zinda * MAX( rzero , SIGN( rone, - phicif(ji,jj) ) ) !subscripts are bad here 127 i1mfr = 1.0 - MAX( rzero , SIGN( rone , - at_i(ji,jj) ) ) 121 128 idfr = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) ) 122 129 iflt = zinda * (1 - i1mfr) * (1 - ifvt ) … … 139 146 140 147 ! computation the solar flux at ocean surface 141 zfcm1(ji,jj) = pfrld(ji,jj) * qsr(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) 148 zfcm1(ji,jj) = pfrld(ji,jj) * qsr(ji,jj) + & 149 & zindb * ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) / MAX( iatte(ji,jj), epsi20 ) 142 150 ! fstric Solar flux transmitted trough the ice 143 151 ! qsr Net short wave heat flux on free ocean 144 152 ! new line 145 fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj)153 fscmbq(ji,jj) = zindb * ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) / MAX( iatte(ji,jj), epsi20 ) 146 154 147 155 ! computation the non solar heat flux at ocean surface … … 178 186 179 187 !!gm this IF prevents the vertorisation of the whole loop 180 IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 181 WRITE(numout,*) ' lim_sbc : heat fluxes ' 182 WRITE(numout,*) ' qsr : ', qsr(jiindx,jjindx) 183 WRITE(numout,*) ' zfcm1 : ', zfcm1(jiindx,jjindx) 184 WRITE(numout,*) ' pfrld : ', pfrld(jiindx,jjindx) 185 WRITE(numout,*) ' fstric : ', fstric (jiindx,jjindx) 186 WRITE(numout,*) 187 WRITE(numout,*) ' qns : ', qns(jiindx,jjindx) 188 WRITE(numout,*) ' zfcm2 : ', zfcm2(jiindx,jjindx) 189 WRITE(numout,*) ' zfcm1 : ', zfcm1(jiindx,jjindx) 190 WRITE(numout,*) ' ifral : ', ifral 191 WRITE(numout,*) ' ial : ', ial 192 WRITE(numout,*) ' qcmif : ', qcmif(jiindx,jjindx) 193 WRITE(numout,*) ' qldif : ', qldif(jiindx,jjindx) 194 WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice 195 WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice 196 WRITE(numout,*) ' ifrdv : ', ifrdv 197 WRITE(numout,*) ' qfvbq : ', qfvbq(jiindx,jjindx) 198 WRITE(numout,*) ' qdtcn : ', qdtcn(jiindx,jjindx) 199 WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice 200 WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice 201 WRITE(numout,*) ' ' 202 WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx) 203 WRITE(numout,*) ' fhmec : ', fhmec(jiindx,jjindx) 204 WRITE(numout,*) ' fheat_rpo : ', fheat_rpo(jiindx,jjindx) 205 WRITE(numout,*) ' fhbri : ', fhbri(jiindx,jjindx) 206 WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 207 ENDIF 188 ! IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 189 ! WRITE(numout,*) 'lim_sbc : heat fluxes ' 190 ! WRITE(numout,*) ' at_i : ', at_i(jiindx,jjindx) 191 ! WRITE(numout,*) ' ht_i : ', SUM( ht_i(jiindx,jjindx,1:jpl) ) 192 ! WRITE(numout,*) ' ht_s : ', SUM( ht_s(jiindx,jjindx,1:jpl) ) 193 ! WRITE(numout,*) 194 ! WRITE(numout,*) ' qsr : ', qsr(jiindx,jjindx) 195 ! WRITE(numout,*) ' zfcm1 : ', zfcm1(jiindx,jjindx) 196 ! WRITE(numout,*) ' pfrld : ', pfrld(jiindx,jjindx) 197 ! WRITE(numout,*) ' fstric : ', fstric (jiindx,jjindx) 198 ! WRITE(numout,*) 199 ! WRITE(numout,*) ' qns : ', qns(jiindx,jjindx) 200 ! WRITE(numout,*) ' zfcm2 : ', zfcm2(jiindx,jjindx) 201 ! WRITE(numout,*) ' zfcm1 : ', zfcm1(jiindx,jjindx) 202 ! WRITE(numout,*) ' ifral : ', ifral 203 ! WRITE(numout,*) ' ial : ', ial 204 ! WRITE(numout,*) ' qcmif : ', qcmif(jiindx,jjindx) 205 ! WRITE(numout,*) ' qldif : ', qldif(jiindx,jjindx) 206 ! !WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice 207 ! !WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice 208 ! WRITE(numout,*) ' ifrdv : ', ifrdv 209 ! WRITE(numout,*) ' qfvbq : ', qfvbq(jiindx,jjindx) 210 ! WRITE(numout,*) ' qdtcn : ', qdtcn(jiindx,jjindx) 211 ! !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice 212 ! !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice 213 ! WRITE(numout,*) ' ' 214 ! WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx) 215 ! WRITE(numout,*) ' fhmec : ', fhmec(jiindx,jjindx) 216 ! WRITE(numout,*) ' fheat_rpo : ', fheat_rpo(jiindx,jjindx) 217 ! WRITE(numout,*) ' fhbri : ', fhbri(jiindx,jjindx) 218 ! WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 219 ! ENDIF 208 220 !!gm end 209 221 END DO … … 236 248 ! computing salt exchanges at the ice/ocean interface 237 249 ! sice should be the same as computed with the ice model 238 zfons = ( soce_0(ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice250 !zfons = ( soce_0(ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice 239 251 ! SOCE 240 zfons = ( sss_m (ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice241 252 !zfons = ( sss_m (ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice 253 zfmm = rdmicif(ji,jj) * r1_rdtice ! IOVINO 242 254 !CT useless ! salt flux for constant salinity 243 255 !CT useless fsalt(ji,jj) = zfons / ( sss_m(ji,jj) + epsi16 ) + fsalt_res(ji,jj) … … 247 259 fsbri(ji,jj) = zinda*fsbri(ji,jj) 248 260 ! converting the salt fluxes from ice to a freshwater flux from ocean 249 fsalt_res(ji,jj) = fsalt_res(ji,jj) / ( sss_m(ji,jj) + epsi16 )250 fseqv(ji,jj) = fseqv(ji,jj) / ( sss_m(ji,jj) + epsi16 )251 fsbri(ji,jj) = fsbri(ji,jj) / ( sss_m(ji,jj) + epsi16 )252 fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) / ( sss_m(ji,jj) + epsi16 )261 ! fsalt_res(ji,jj) = fsalt_res(ji,jj) / ( sss_m(ji,jj) + epsi16 ) 262 ! fseqv(ji,jj) = fseqv(ji,jj) / ( sss_m(ji,jj) + epsi16 ) 263 ! fsbri(ji,jj) = fsbri(ji,jj) / ( sss_m(ji,jj) + epsi16 ) 264 ! fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) / ( sss_m(ji,jj) + epsi16 ) 253 265 254 266 ! freshwater mass exchange (positive to the ice, negative for the ocean ?) … … 258 270 ! POSITIVE FRESHWATER FLUX FROM THE OCEAN TO THE ICE [kg.m-2.s-1] 259 271 260 emp(ji,jj) = - zpme 272 emp(ji,jj) = - zpme + zfmm ! volume flux IOVINO 273 ! emp(ji,jj) = - zpme 261 274 END DO 262 275 END DO 263 276 264 277 IF( num_sal == 2 ) THEN ! variable ice salinity: brine drainage included in the salt flux 265 emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)278 emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ! + emp(:,:) ! IOVINO 266 279 ELSE ! constant ice salinity: 267 emps(:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)280 emps(:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ! + emp(:,:) ! IOVINO 268 281 ENDIF 282 269 283 270 284 !-----------------------------------------------! … … 402 416 END WHERE 403 417 ENDIF 418 ! clem modif 419 iatte(:,:) = 1._wp 420 oatte(:,:) = 1._wp 404 421 ! 405 422 END SUBROUTINE lim_sbc_init -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r3294 r3938 11 11 !! 3.3 ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas) 12 12 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 13 !! - ! 2012-05 (C. Rousset) add penetration solar flux 13 14 !!---------------------------------------------------------------------- 14 15 #if defined key_lim3 … … 39 40 USE in_out_manager ! I/O manager 40 41 USE prtctl ! Print control 42 USE lib_fortran ! to use key_nosignedzero 41 43 42 44 IMPLICIT NONE … … 91 93 REAL(wp) :: zfntlat, zpareff, zareamin, zcoef ! - - 92 94 REAL(wp), POINTER, DIMENSION(:,:) :: zqlbsbq ! link with lead energy budget qldif 95 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 93 96 !!------------------------------------------------------------------- 94 97 95 98 CALL wrk_alloc( jpi, jpj, zqlbsbq ) 96 99 100 ! ------------------------------- 101 !- check conservation (C Rousset) 102 IF (ln_limdiahsb) THEN 103 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 104 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 105 zchk_fw_b = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) 106 zchk_fs_b = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) 107 ENDIF 108 !- check conservation (C Rousset) 109 ! ------------------------------- 110 97 111 !------------------------------------------------------------------------------! 98 112 ! 1) Initialization of diagnostic variables ! … … 108 122 DO ji = 1, jpi 109 123 !Energy of melting q(S,T) [J.m-3] 110 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * nlay_i124 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * REAL( nlay_i ) 111 125 !0 if no ice and 1 if yes 112 126 zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i(ji,jj,jl) ) ) … … 120 134 DO ji = 1, jpi 121 135 !Energy of melting q(S,T) [J.m-3] 122 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * nlay_s136 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL( nlay_s ) 123 137 !0 if no ice and 1 if yes 124 138 zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s(ji,jj,jl) ) ) … … 133 147 ! 1.3) Set some dummies to 0 134 148 !----------------------------- 135 rdvosif(:,:) = 0.e0 ! variation of ice volume at surface136 rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom137 fdvolif(:,:) = 0.e0 ! total variation of ice volume138 rdvonif(:,:) = 0.e0 ! lateral variation of ice volume139 fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice140 ffltbif(:,:) = 0.e0 ! linked with fstric141 qfvbq (:,:) = 0.e0 ! linked with fstric142 rdmsnif(:,:) = 0.e0 ! variation of snow mass per unit area143 rdmicif(:,:) = 0.e0 ! variation of ice mass per unit area144 hicifp (:,:) = 0.e0 ! daily thermodynamic ice production.145 fsbri (:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean146 fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean147 fseqv (:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay149 !clem rdvosif(:,:) = 0.e0 ! variation of ice volume at surface 150 !clem rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom 151 !clem fdvolif(:,:) = 0.e0 ! total variation of ice volume 152 !clem rdvonif(:,:) = 0.e0 ! lateral variation of ice volume 153 !clem fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice 154 !clem ffltbif(:,:) = 0.e0 ! linked with fstric 155 !clem qfvbq (:,:) = 0.e0 ! linked with fstric 156 !clem rdmsnif(:,:) = 0.e0 ! variation of snow mass per unit area 157 !clem rdmicif(:,:) = 0.e0 ! variation of ice mass per unit area 158 !clem hicifp (:,:) = 0.e0 ! daily thermodynamic ice production. 159 !clem fsbri (:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean 160 !clem fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean 161 !clem fseqv (:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay 148 162 149 163 !----------------------------------- … … 164 178 !CDIR NOVERRCHK 165 179 DO ji = 1, jpi 166 zthsnice = SUM( ht_s(ji,jj,1:jpl) ) + SUM( ht_i(ji,jj,1:jpl) )167 zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice) ) )180 !clem zthsnice = SUM( ht_s(ji,jj,1:jpl) ) + SUM( ht_i(ji,jj,1:jpl) ) 181 !clem zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice + epsi20 ) ) ) 168 182 phicif(ji,jj) = vt_i(ji,jj) 169 183 pfrld(ji,jj) = 1.0 - at_i(ji,jj) 170 zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )184 zinda = tms(ji,jj) * (1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) ) 171 185 ! 172 186 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 179 193 180 194 ! here the drag will depend on ice thickness and type (0.006) 181 fdtcn(ji,jj) = zind b* rau0 * rcp * 0.006 * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) )195 fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) ) 182 196 ! also category dependent 183 197 ! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead 184 qdtcn(ji,jj) = zind b* fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice198 qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 185 199 ! 186 200 ! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) 187 201 ! ! caution: exponent betas used as more snow can fallinto leads 188 202 qldif(ji,jj) = tms(ji,jj) * rdt_ice * ( & 189 & pfrld(ji,jj) * ( qsr(ji,jj) & ! solar heat203 & pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif 190 204 & + qns(ji,jj) & ! non solar heat 191 205 & + fdtcn(ji,jj) & ! turbulent ice-ocean heat 192 & + fsbbq(ji,jj) * ( 1.0 - zind b) ) & ! residual heat from previous step206 & + fsbbq(ji,jj) * ( 1.0 - zinda ) ) & ! residual heat from previous step 193 207 & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus ) ! latent heat of sprecip melting 194 208 ! … … 205 219 ! 206 220 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 207 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda )221 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) !!clem * ( 1. - zinda ) 208 222 ! 209 223 ! oceanic heat flux (limthd_dh) 210 fbif (ji,jj) = zind b* ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) )224 fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 211 225 ! 212 226 END DO … … 294 308 CALL tab_2d_1d( nbpb, fstbif_1d (1:nbpb), fstric , jpi, jpj, npb(1:nbpb) ) 295 309 CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb), qfvbq , jpi, jpj, npb(1:nbpb) ) 296 310 CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte(:,:) , jpi, jpj, npb(1:nbpb) ) ! clem modif 311 CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte(:,:) , jpi, jpj, npb(1:nbpb) ) ! clem modif 297 312 !-------------------------------- 298 313 ! 4.3) Thermodynamic processes … … 305 320 CALL lim_thd_dif( 1, nbpb, jl ) ! Ice/Snow Temperature profile ! 306 321 ! !---------------------------------! 307 322 ! 308 323 CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting compulsory for limthd_dh 309 324 ! 310 325 IF( con_i ) CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 311 326 IF( con_i ) CALL lim_thd_con_dif( 1 , nbpb , jl ) … … 314 329 CALL lim_thd_dh( 1, nbpb, jl ) ! Ice/Snow thickness ! 315 330 ! !---------------------------------! 316 317 331 ! !---------------------------------! 318 332 CALL lim_thd_ent( 1, nbpb, jl ) ! Ice/Snow enthalpy remapping ! 319 333 ! !---------------------------------! 320 321 334 ! !---------------------------------! 322 335 CALL lim_thd_sal( 1, nbpb ) ! Ice salinity computation ! 323 336 ! !---------------------------------! 324 325 337 ! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 326 338 IF( con_i ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) … … 418 430 ! 5.4) Diagnostic thermodynamic growth rates 419 431 !-------------------------------------------- 420 d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes421 432 !clem@useless d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes 433 !clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) / rdt_ice * 86400.0 422 434 423 435 IF( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) … … 455 467 ENDIF 456 468 ! 469 ! ------------------------------- 470 !- check conservation (C Rousset) 471 IF (ln_limdiahsb) THEN 472 zchk_fs = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 473 zchk_fw = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 474 475 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 476 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 477 478 IF(lwp) THEN 479 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limthd) = ',(zchk_v_i * 86400.) 480 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * 86400.) 481 IF ( MINVAL( v_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limthd) = ',(MINVAL(v_i) * 1.e-3) 482 IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limthd) = ',MAXVAL(SUM(a_i,dim=3)) 483 ENDIF 484 ENDIF 485 !- check conservation (C Rousset) 486 ! ------------------------------- 487 ! 457 488 CALL wrk_dealloc( jpi, jpj, zqlbsbq ) 458 489 ! … … 479 510 DO jk = 1, nlay_i ! total q over all layers, ice [J.m-2] 480 511 DO ji = kideb, kiut 481 etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i512 etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 482 513 eti (ji,jl) = eti(ji,jl) + etilayer(ji,jk) 483 514 END DO 484 515 END DO 485 516 DO ji = kideb, kiut ! total q over all layers, snow [J.m-2] 486 ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / nlay_s517 ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / REAL( nlay_s ) 487 518 END DO 488 519 ! … … 799 830 !!------------------------------------------------------------------- 800 831 NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, & 801 & hicmin, hiclim, amax ,&832 & hicmin, hiclim, & 802 833 & sbeta , parlat, hakspl, hibspl, exld, & 803 834 & hakdif, hnzst , thth , parsub, alphs, betas, & … … 825 856 WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin 826 857 WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim 827 WRITE(numout,*)' maximum lead fraction amax = ', amax828 858 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 829 859 WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r3294 r3938 24 24 USE lib_mpp ! MPP library 25 25 USE wrk_nemo ! work arrays 26 USE lib_fortran ! to use key_nosignedzero 26 27 27 28 IMPLICIT NONE … … 72 73 INTEGER :: ji , jk ! dummy loop indices 73 74 INTEGER :: zji, zjj ! 2D corresponding indices to ji 74 INTEGER :: isnow ! switch for presence (1) or absence (0) of snow75 INTEGER :: isnowic ! snow ice formation not76 INTEGER :: i_ice_switch ! ice thickness above a certain treshold or not77 75 INTEGER :: iter 76 77 REAL(wp) :: isnow ! switch for presence (1) or absence (0) of snow 78 REAL(wp) :: isnowic ! snow ice formation not 79 REAL(wp) :: i_ice_switch ! ice thickness above a certain treshold or not 78 80 79 81 REAL(wp) :: zzfmass_i, zihgnew ! local scalar … … 118 120 REAL(wp), POINTER, DIMENSION(:,:) :: zqt_i_lay ! total ice heat content 119 121 122 ! mass and salt flux (clem) 123 REAL(wp) :: zdvres 124 REAL(wp), POINTER, DIMENSION(:) :: zviold, zvsold ! old ice volume... 125 120 126 ! Heat conservation 121 127 INTEGER :: num_iter_max, numce_dh … … 130 136 CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) 131 137 138 CALL wrk_alloc( jpij, zviold, zvsold ) ! clem 139 132 140 zfsalt_melt(:) = 0._wp 133 141 ftotal_fin(:) = 0._wp 134 142 zfdt_init(:) = 0._wp 135 143 zfdt_final(:) = 0._wp 144 dh_i_surf(:) = 0._wp 145 dh_i_bott(:) = 0._wp 146 dh_snowice(:) = 0._wp 136 147 137 148 DO ji = kideb, kiut 138 149 old_ht_i_b(ji) = ht_i_b(ji) 139 150 old_ht_s_b(ji) = ht_s_b(ji) 151 zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem 152 zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem 140 153 END DO 141 154 ! … … 145 158 ! 146 159 DO ji = kideb, kiut 147 isnow = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji)) ) )160 isnow = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) 148 161 ztfs(ji) = isnow * rtt + ( 1.0 - isnow ) * rtt 149 162 z_f_surf(ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) … … 162 175 ! 163 176 DO ji = kideb, kiut ! Layer thickness 164 zh_i(ji) = ht_i_b(ji) / nlay_i165 zh_s(ji) = ht_s_b(ji) / nlay_s177 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 178 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 166 179 END DO 167 180 ! … … 169 182 DO jk = 1, nlay_s 170 183 DO ji = kideb, kiut 171 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s184 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) 172 185 END DO 173 186 END DO … … 176 189 DO jk = 1, nlay_i 177 190 DO ji = kideb, kiut 178 zzc = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i191 zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 179 192 zqt_i(ji) = zqt_i(ji) + zzc 180 193 zqt_i_lay(ji,jk) = zzc … … 243 256 ht_s_b(ji) = MAX( zzero , zhsnew ) 244 257 ! Volume and mass variations of snow 245 dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) ) 258 ! dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) ) 259 dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) ! IOVINO 246 260 dvsbq_1d (ji) = MIN( zzero, dvsbq_1d(ji) ) 247 rdmsnif_1d(ji) = rdmsnif_1d(ji) + rhosn * dvsbq_1d(ji)261 !clem rdmsnif_1d(ji) = rdmsnif_1d(ji) + rhosn * dvsbq_1d(ji) 248 262 END DO ! ji 249 263 … … 252 266 !-------------------------- 253 267 DO ji = kideb, kiut 254 dh_i_surf(ji) = 0._wp255 268 z_f_surf (ji) = zqfont_su(ji) / rdt_ice ! heat conservation test 256 269 zdq_i (ji) = 0._wp … … 270 283 zdq_i (ji) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 271 284 ! 272 ! contribution to ice-ocean salt flux273 zji = MOD( npb(ji) - 1 , jpi ) + 1274 zjj = ( npb(ji) - 1 ) / jpi + 1275 zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji)) * a_i_b(ji) &276 & * MIN( zdeltah(ji,jk) , 0. e0 ) * rhoic / rdt_ice285 ! IOVINO 286 !zfsalt_melt(ji) = zfsalt_melt(ji) - sm_i_b(ji) * a_i_b(ji) & 287 ! & * MIN( zdeltah(ji,jk) , 0.e0 ) * rhoic / rdt_ice 288 fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 289 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 277 290 END DO 278 291 END DO … … 331 344 DO ji = kideb,kiut 332 345 q_s_b (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 333 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s! heat conservation346 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) ! heat conservation 334 347 END DO 335 348 END DO … … 372 385 ! Basal growth rate = - F*dt / q 373 386 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 387 fseqv_1d(ji) = fseqv_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic / rdt_ice 374 388 ENDIF 375 389 END DO … … 437 451 dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) ) & 438 452 & / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji) 453 fseqv_1d(ji) = fseqv_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic / rdt_ice 439 454 ENDIF ! heat budget 440 455 END DO … … 473 488 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) 474 489 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 475 ! contribution to salt flux476 zji = MOD( npb(ji) - 1, jpi ) + 1477 zjj = ( npb(ji) - 1 ) / jpi + 1478 zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji) ) * a_i_b(ji) &479 & * MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice480 490 ENDIF 491 ! IOVINO contribution to salt flux 492 !zfsalt_melt(ji) = zfsalt_melt(ji) - sm_i_b(ji) * a_i_b(ji) & 493 ! & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 494 fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 495 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 481 496 ENDIF 482 497 END DO ! ji … … 529 544 ELSE ; zdhbf = dh_i_bott(ji) 530 545 ENDIF 546 zdvres = zdhbf - dh_i_bott(ji) 547 dh_i_bott(ji) = zdhbf 548 fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic / rdt_ice 531 549 ! ! excessive energy is sent to lateral ablation 532 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) & 533 & * ( zdhbf - dh_i_bott(ji) ) / rdt_ice 534 dh_i_bott(ji) = zdhbf 550 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres / rdt_ice 535 551 ! !since ice volume is only used for outputs, we keep it global for all categories 536 552 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 537 553 ! !new ice thickness 538 554 zhgnew (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 539 ! ! diagnostic ( bottom ice growth )540 zji = MOD( npb(ji) - 1, jpi ) + 1541 zjj = ( npb(ji) - 1 ) / jpi + 1542 diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice543 diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) / rdt_ice544 diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice545 555 END DO 546 556 … … 554 564 ! Adapt the remaining energy if too much ice melts 555 565 !-------------------------------------------------- 566 zdvres = MAX( 0._wp, - zhgnew(ji) ) 556 567 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 557 ! 0 if no more ice568 dh_i_bott (ji) = dh_i_bott(ji) + zdvres ! clem@bug 558 569 zhgnew (ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 559 570 ! remaining heat … … 581 592 ! 582 593 ! ! mass variation cumulated over category 583 rdmsnif_1d(ji) = rdmsnif_1d(ji) + zzfmass_s ! snow584 rdmicif_1d(ji) = rdmicif_1d(ji) + zzfmass_i ! ice594 !clem rdmsnif_1d(ji) = rdmsnif_1d(ji) + zzfmass_s ! snow 595 !clem rdmicif_1d(ji) = rdmicif_1d(ji) + zzfmass_i ! ice 585 596 586 597 ! Remaining heat to the ocean 587 598 !--------------------------------- 588 599 focea(ji) = - zfdt_final(ji) / rdt_ice ! focea is in W.m-2 * dt 589 600 ! salt flux 601 fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic / rdt_ice 602 ! ! diagnostic ( bottom ice growth ) 603 zji = MOD( npb(ji) - 1, jpi ) + 1 604 zjj = ( npb(ji) - 1 ) / jpi + 1 605 diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice 606 diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) / rdt_ice 607 diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice 590 608 END DO 591 609 … … 602 620 zjj = ( npb(ji) - 1 ) / jpi + 1 603 621 ! new lines 604 IF( num_sal == 4 ) THEN 605 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) & 606 & + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice 607 ELSE 608 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) & 609 & + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice 610 ENDIF 622 !IF( num_sal == 4 ) THEN 623 ! ! IOVINO 624 ! fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) & 625 ! & - (1.0 - zihgnew) * zfmass_i(ji) * bulk_sal / rdt_ice 626 !ELSE 627 ! ! IOVINO 628 ! fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) & 629 ! & - (1.0 - zihgnew) * zfmass_i(ji) * sm_i_b(ji) / rdt_ice 630 !ENDIF 611 631 ! Heat flux 612 632 ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1 … … 619 639 qldif_1d(ji) = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice & 620 640 & + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice 641 !IF ( (zji.eq.jiindx).AND.(zjj.eq.jjindx) ) THEN 642 ! !clemclem 643 ! WRITE(numout,*) 'lim_thd_dh : qldif = ', qldif_1d(ji) 644 ! !clemclem 645 !ENDIF 621 646 END DO ! ji 622 647 … … 656 681 dmgwi_1d (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 657 682 658 rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic659 rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn683 !clem rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 684 !clem rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 660 685 661 686 ! Equivalent salt flux (1) Snow-ice formation component … … 667 692 ELSE ; zsm_snowice = ( rhoic - rhosn ) / rhoic * sss_m(zji,zjj) 668 693 ENDIF 669 IF( num_sal == 4 ) THEN670 fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal ) * a_i_b(ji) &671 & * ( zhgnew(ji) - ht_i_b(ji)) * rhoic / rdt_ice672 ELSE673 fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - zsm_snowice ) * a_i_b(ji) &674 & * ( zhgnew(ji) - ht_i_b(ji)) * rhoic / rdt_ice675 ENDIF694 !IF( num_sal == 4 ) THEN 695 ! ! IOVINO 696 ! fseqv_1d(ji) = fseqv_1d(ji) - bulk_sal * a_i_b(ji) * dh_snowice(ji) * rhoic / rdt_ice 697 !ELSE 698 ! ! IOVINO 699 ! fseqv_1d(ji) = fseqv_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic / rdt_ice 700 !ENDIF 676 701 ! entrapment during snow ice formation 677 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 678 isnowic = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch 679 IF( num_sal == 2 .OR. num_sal == 4 ) & 680 dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) & 681 & + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13) & 682 & - sm_i_b(ji) ) * isnowic 702 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + epsi13 ) ) 703 isnowic = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch 704 705 !clem IF( num_sal == 2 .OR. num_sal == 4 ) & 706 !clem dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) & 707 !clem & + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13) & 708 !clem & - sm_i_b(ji) ) * isnowic 709 IF ( num_sal == 2 .OR. num_sal == 4 ) & 710 & dsm_i_si_1d(ji) = ( ( zsm_snowice * dh_snowice(ji) + sm_i_b(ji) * ht_i_b(ji) ) & 711 & / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 ) - sm_i_b(ji) ) * isnowic 683 712 684 713 ! Actualize new snow and ice thickness. … … 693 722 zjj = ( npb(ji) - 1 ) / jpi + 1 694 723 diag_sni_gr(zji,zjj) = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / rdt_ice 695 ! 724 725 ! salt flux 726 fseqv_1d(ji) = fseqv_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic / rdt_ice 727 !-------------------------------- 728 ! Update mass fluxes (clem) 729 !-------------------------------- 730 rdmicif_1d(ji) = rdmicif_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic 731 rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn 732 696 733 END DO !ji 697 734 ! … … 700 737 CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) 701 738 CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) 739 ! 740 CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem 702 741 ! 703 742 END SUBROUTINE lim_thd_dh -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r3351 r3938 10 10 !! ! 04-2007 (M. Vancoppenolle) Energy conservation 11 11 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 12 !! - ! 2012-05 (C. Rousset) add penetration solar flux 12 13 !!---------------------------------------------------------------------- 13 14 #if defined key_lim3 … … 23 24 USE lib_mpp ! MPP library 24 25 USE wrk_nemo ! work arrays 26 USE lib_fortran ! to use key_nosignedzero 25 27 26 28 IMPLICIT NONE … … 102 104 INTEGER :: nconv ! number of iterations in iterative procedure 103 105 INTEGER :: minnumeqmin, maxnumeqmax 104 INTEGER, DIMENSION(kiut) :: numeqmin ! reference number of top equation 105 INTEGER, DIMENSION(kiut) :: numeqmax ! reference number of bottom equation 106 INTEGER, DIMENSION(kiut) :: isnow ! switch for presence (1) or absence (0) of snow 106 107 INTEGER , DIMENSION(kiut) :: numeqmin ! reference number of top equation 108 INTEGER , DIMENSION(kiut) :: numeqmax ! reference number of bottom equation 109 INTEGER , DIMENSION(kiut) :: isnow ! switch for presence (1) or absence (0) of snow 110 111 !! * New local variables 112 REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztcond_i !Ice thermal conductivity 113 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zradtr_i !Radiation transmitted through the ice 114 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zradab_i !Radiation absorbed in the ice 115 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zkappa_i !Kappa factor in the ice 116 117 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zradtr_s !Radiation transmited through the snow 118 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zradab_s !Radiation absorbed in the snow 119 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zkappa_s !Kappa factor in the snow 120 121 REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztiold !Old temperature in the ice 122 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zeta_i !Eta factor in the ice 123 REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztitemp !Temporary temperature in the ice to check the convergence 124 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zspeche_i !Ice specific heat 125 REAL(wp), DIMENSION(kiut,0:nlay_i) :: z_i !Vertical cotes of the layers in the ice 126 127 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zeta_s !Eta factor in the snow 128 REAL(wp), DIMENSION(kiut,0:nlay_s) :: ztstemp !Temporary temperature in the snow to check the convergence 129 REAL(wp), DIMENSION(kiut,0:nlay_s) :: ztsold !Temporary temperature in the snow 130 REAL(wp), DIMENSION(kiut,0:nlay_s) :: z_s !Vertical cotes of the layers in the snow 131 132 REAL(wp), DIMENSION(kiut,jkmax+2) :: zindterm ! Independent term 133 REAL(wp), DIMENSION(kiut,jkmax+2) :: zindtbis ! temporary independent term 134 REAL(wp), DIMENSION(kiut,jkmax+2) :: zdiagbis 135 REAL(wp), DIMENSION(kiut,jkmax+2,3) :: ztrid ! tridiagonal system terms 136 137 REAL(wp), DIMENSION(kiut) :: ztfs ! ice melting point 138 REAL(wp), DIMENSION(kiut) :: ztsuold ! old surface temperature (before the iterative procedure ) 139 REAL(wp), DIMENSION(kiut) :: ztsuoldit ! surface temperature at previous iteration 140 REAL(wp), DIMENSION(kiut) :: zh_i ! ice layer thickness 141 REAL(wp), DIMENSION(kiut) :: zh_s ! snow layer thickness 142 REAL(wp), DIMENSION(kiut) :: zfsw ! solar radiation absorbed at the surface 143 REAL(wp), DIMENSION(kiut) :: zf ! surface flux function 144 REAL(wp), DIMENSION(kiut) :: dzf ! derivative of the surface flux function 145 107 146 REAL(wp) :: zeps = 1.e-10_wp ! 108 147 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system … … 113 152 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 114 153 REAL(wp) :: zht_smin = 1.e-4_wp ! minimum snow depth 154 115 155 REAL(wp) :: ztmelt_i ! ice melting temperature 116 156 REAL(wp) :: zerritmax ! current maximal error on temperature 117 REAL(wp), DIMENSION(kiut) :: ztfs ! ice melting point 118 REAL(wp), DIMENSION(kiut) :: ztsuold ! old surface temperature (before the iterative procedure ) 119 REAL(wp), DIMENSION(kiut) :: ztsuoldit ! surface temperature at previous iteration 120 REAL(wp), DIMENSION(kiut) :: zh_i ! ice layer thickness 121 REAL(wp), DIMENSION(kiut) :: zh_s ! snow layer thickness 122 REAL(wp), DIMENSION(kiut) :: zfsw ! solar radiation absorbed at the surface 123 REAL(wp), DIMENSION(kiut) :: zf ! surface flux function 124 REAL(wp), DIMENSION(kiut) :: dzf ! derivative of the surface flux function 125 REAL(wp), DIMENSION(kiut) :: zerrit ! current error on temperature 126 REAL(wp), DIMENSION(kiut) :: zdifcase ! case of the equation resolution (1->4) 127 REAL(wp), DIMENSION(kiut) :: zftrice ! solar radiation transmitted through the ice 157 REAL(wp), DIMENSION(kiut) :: zerrit ! current error on temperature 158 REAL(wp), DIMENSION(kiut) :: zdifcase ! case of the equation resolution (1->4) 159 REAL(wp), DIMENSION(kiut) :: zftrice ! solar radiation transmitted through the ice 128 160 REAL(wp), DIMENSION(kiut) :: zihic, zhsu 129 REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztcond_i ! Ice thermal conductivity130 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice131 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice132 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zkappa_i ! Kappa factor in the ice133 REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztiold ! Old temperature in the ice134 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zeta_i ! Eta factor in the ice135 REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztitemp ! Temporary temperature in the ice to check the convergence136 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zspeche_i ! Ice specific heat137 REAL(wp), DIMENSION(kiut,0:nlay_i) :: z_i ! Vertical cotes of the layers in the ice138 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow139 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow140 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zkappa_s ! Kappa factor in the snow141 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zeta_s ! Eta factor in the snow142 REAL(wp), DIMENSION(kiut,0:nlay_s) :: ztstemp ! Temporary temperature in the snow to check the convergence143 REAL(wp), DIMENSION(kiut,0:nlay_s) :: ztsold ! Temporary temperature in the snow144 REAL(wp), DIMENSION(kiut,0:nlay_s) :: z_s ! Vertical cotes of the layers in the snow145 REAL(wp), DIMENSION(kiut,jkmax+2) :: zindterm ! Independent term146 REAL(wp), DIMENSION(kiut,jkmax+2) :: zindtbis ! temporary independent term147 REAL(wp), DIMENSION(kiut,jkmax+2) :: zdiagbis148 REAL(wp), DIMENSION(kiut,jkmax+2,3) :: ztrid ! tridiagonal system terms149 161 !!------------------------------------------------------------------ 150 151 ! 162 ! 152 163 !------------------------------------------------------------------------------! 153 164 ! 1) Initialization ! … … 156 167 DO ji = kideb , kiut 157 168 ! is there snow or not 158 isnow(ji)= INT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) )169 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) ) 159 170 ! surface temperature of fusion 160 171 !!gm ??? ztfs(ji) = rtt !!!???? 161 ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt172 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 162 173 ! layer thickness 163 zh_i(ji) = ht_i_b(ji) / nlay_i164 zh_s(ji) = ht_s_b(ji) / nlay_s174 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 175 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 165 176 END DO 166 177 … … 174 185 DO layer = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 175 186 DO ji = kideb , kiut 176 z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s187 z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s ) 177 188 END DO 178 189 END DO … … 180 191 DO layer = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 181 192 DO ji = kideb , kiut 182 z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i193 z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i ) 183 194 END DO 184 195 END DO … … 201 212 DO ji = kideb , kiut 202 213 ! switches 203 isnow(ji) = INT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ) )214 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ) ) 204 215 ! hs > 0, isnow = 1 205 216 zhsu (ji) = hnzst ! threshold for the computation of i0 206 217 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) ) 207 218 208 i0(ji) = ( 1._wp- isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )219 i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 209 220 !fr1_i0_1d = i0 for a thin ice surface 210 221 !fr1_i0_2d = i0 for a thick ice surface … … 243 254 244 255 DO ji = kideb, kiut ! ice initialization 245 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp- isnow(ji) )256 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1 - isnow(ji) ) 246 257 END DO 247 258 … … 256 267 257 268 DO ji = kideb, kiut ! Radiation transmitted below the ice 258 fstbif_1d(ji) = fstbif_1d(ji) + zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji)269 fstbif_1d(ji) = fstbif_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif 259 270 END DO 260 271 … … 264 275 ii = MOD( npb(ji) - 1, jpi ) + 1 265 276 ij = ( npb(ji) - 1 ) / jpi + 1 266 fstroc(ii,ij,jl) = zradtr_i(ji,nlay_i)277 fstroc(ii,ij,jl) = iatte_1d(ji) * zradtr_i(ji,nlay_i) ! clem modif 267 278 END DO 268 279 ! +++++ … … 273 284 END DO 274 285 END DO 275 276 286 277 287 ! … … 377 387 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & 378 388 (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 379 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)* isnow(ji) &380 + zkappa_i(ji,0)* (1.0-isnow(ji))389 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 390 + zkappa_i(ji,0)*REAL( 1 - isnow(ji) ) 381 391 END DO 382 392 ! … … 659 669 t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 660 670 * t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) & 661 * MAX(0.0,SIGN(1.0,ht_s_b(ji) -zeps))671 * MAX(0.0,SIGN(1.0,ht_s_b(ji))) 662 672 663 673 ! surface temperature 664 isnow(ji) = INT(1.0-max(0.0,sign(1.0,-ht_s_b(ji))))674 isnow(ji) = NINT(1.0-max(0.0,SIGN(1.0,-ht_s_b(ji)))) 665 675 ztsuoldit(ji) = t_su_b(ji) 666 676 IF (t_su_b(ji) .LT. ztfs(ji)) & 667 t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( isnow(ji)*t_s_b(ji,1) &668 & + (1.0-isnow(ji))*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))677 t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) & 678 & + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 669 679 END DO 670 680 ! … … 718 728 DO ji = kideb, kiut 719 729 ! ! update of latent heat fluxes 720 qla_ice_1d (ji) = qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) 730 ! qla_ice_1d (ji) = qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji)! - ztsuold(ji) ) ! TECLIM change 731 qla_ice_1d (ji) = MAX( 0.e0, qla_ice_1d (ji) + dqla_ice_1d(ji) * (t_su_b(ji) - ztsuold(ji) ) ) 721 732 ! ! surface ice conduction flux 722 isnow(ji) = INT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) )723 fc_su(ji) = - isnow(ji)* zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) &724 & - ( 1._wp- isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji))733 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) ) 734 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) & 735 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji)) 725 736 ! ! bottom ice conduction flux 726 737 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) … … 733 744 DO ji = kideb, kiut 734 745 ! Upper snow value 735 fc_s(ji,0) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) )746 fc_s(ji,0) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) ) 736 747 ! Bott. snow value 737 fc_s(ji,1) = - isnow(ji)* zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) )748 fc_s(ji,1) = - REAL( isnow(ji) ) * zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) ) 738 749 END DO 739 750 DO ji = kideb, kiut ! Upper ice layer 740 fc_i(ji,0) = - isnow(ji) * & ! interface flux if there is snow751 fc_i(ji,0) = - REAL( isnow(ji) ) * & ! interface flux if there is snow 741 752 ( zkappa_i(ji,0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) & 742 - ( 1.0- isnow(ji) ) * ( zkappa_i(ji,0) * &753 - REAL( 1 - isnow(ji) ) * ( zkappa_i(ji,0) * & 743 754 zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not 744 755 END DO … … 755 766 ENDIF 756 767 ! 768 757 769 END SUBROUTINE lim_thd_dif 758 770 -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r3294 r3938 29 29 USE lib_mpp ! MPP library 30 30 USE wrk_nemo ! work arrays 31 USE lib_fortran ! to use key_nosignedzero 31 32 32 33 IMPLICIT NONE … … 144 145 145 146 DO ji = kideb, kiut 146 zh_i(ji) = old_ht_i_b(ji) / nlay_i147 zh_s(ji) = old_ht_s_b(ji) / nlay_s147 zh_i(ji) = old_ht_i_b(ji) / REAL( nlay_i ) 148 zh_s(ji) = old_ht_s_b(ji) / REAL( nlay_s ) 148 149 END DO 149 150 … … 165 166 DO jk = 1, nlays0 166 167 DO ji = kideb, kiut 167 snind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))) &168 + snind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))))168 snind(ji) = jk * NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)))) & 169 + snind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji))))) 169 170 zdeltah(ji)= zdeltah(ji) + zh_s(ji) 170 171 END DO ! ji … … 174 175 ! 0 if not 175 176 DO ji = kideb, kiut 176 snswi(ji) = MAX(0, INT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji)))))177 snswi(ji) = MAX(0,NINT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji))))) 177 178 END DO ! ji 178 179 … … 189 190 DO jk = 1, nlayi0 190 191 DO ji = kideb, kiut 191 icsuind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))) &192 + icsuind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))))192 icsuind(ji) = jk * NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)))) & 193 + icsuind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji))))) 193 194 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 194 195 END DO ! ji … … 199 200 ! 0 if not 200 201 DO ji = kideb, kiut 201 icsuswi(ji) = MAX(0, INT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) )202 icsuswi(ji) = MAX(0,NINT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) ) 202 203 ENDDO 203 204 … … 215 216 DO jk = nlayi0, 1, -1 216 217 DO ji = kideb, kiut 217 icboind(ji) = (nlayi0+1-jk) * INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))) &218 & + icboind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))))218 icboind(ji) = (nlayi0+1-jk) * NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)))) & 219 & + icboind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji))))) 219 220 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 220 221 END DO … … 231 232 ! 0 if ablation is on the way 232 233 DO ji = kideb, kiut 233 icboswi(ji) = MAX(0, INT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji)))))234 icboswi(ji) = MAX(0,NINT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji))))) 234 235 END DO 235 236 … … 247 248 DO ji = kideb, kiut 248 249 snicind(ji) = (nlays0+1-jk) & 249 * INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))) + snicind(ji) &250 * (1 - INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))))250 * NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)))) + snicind(ji) & 251 * (1 - NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji))))) 251 252 zdeltah(ji) = zdeltah(ji) + zh_s(ji) 252 253 END DO … … 257 258 ! 0 if not 258 259 DO ji = kideb, kiut 259 snicswi(ji) = MAX(0, INT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji)))))260 snicswi(ji) = MAX(0,NINT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji))))) 260 261 ENDDO 261 262 … … 278 279 279 280 DO ji = kideb, kiut 280 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1 .- snicind(ji) ) * snicswi(ji)281 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1 - snicind(ji) ) * snicswi(ji) 281 282 ! cotes of the top of the layers 282 283 zm0(ji,0) = 0._wp … … 290 291 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 291 292 limsum = MIN( limsum , nlay_s ) 292 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * limsum293 END DO 294 END DO 295 296 DO ji = kideb, kiut 297 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - snicswi(ji) * dh_snowice(ji) + zh_s(ji) * nlays0298 zm0(ji,1) = dh_s_tot(ji) * (1 -snswi(ji) ) + snswi(ji) * zm0(ji,1)293 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * REAL( limsum ) 294 END DO 295 END DO 296 297 DO ji = kideb, kiut 298 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - REAL( snicswi(ji) ) * dh_snowice(ji) + zh_s(ji) * REAL( nlays0 ) 299 zm0(ji,1) = dh_s_tot(ji) * REAL( 1 - snswi(ji) ) + REAL( snswi(ji) ) * zm0(ji,1) 299 300 END DO 300 301 … … 308 309 309 310 DO ji = kideb, kiut ! layer heat content 310 qm0 (ji,1) = rhosn * ( cpic * ( rtt - ( 1.- snswi(ji) ) * tatm_ice_1d(ji) &311 & - snswi(ji)* t_s_b (ji,1) ) &311 qm0 (ji,1) = rhosn * ( cpic * ( rtt - REAL( 1 - snswi(ji) ) * tatm_ice_1d(ji) & 312 & - REAL( snswi(ji) ) * t_s_b (ji,1) ) & 312 313 & + lfus ) * zthick0(ji,1) 313 314 zqts_in(ji) = zqts_in(ji) + qm0(ji,1) … … 319 320 limsum = MIN( limsum , nlay_s ) 320 321 qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk) 321 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, epsi20- ht_s_b(ji) ) )322 zqts_in(ji) = zqts_in(ji) + ( 1.- snswi(ji) ) * qm0(ji,jk) * zswitch322 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, - ht_s_b(ji) ) ) 323 zqts_in(ji) = zqts_in(ji) + REAL( 1 - snswi(ji) ) * qm0(ji,jk) * zswitch 323 324 END DO ! jk 324 325 END DO ! ji … … 359 360 !------------------- 360 361 DO ji = kideb, kiut 361 zh_s(ji) = ht_s_b(ji) / nlay_s362 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 362 363 z_s(ji,0) = 0._wp 363 364 ENDDO … … 365 366 DO jk = 1, nlay_s 366 367 DO ji = kideb, kiut 367 z_s(ji,jk) = zh_s(ji) * jk368 z_s(ji,jk) = zh_s(ji) * REAL( jk ) 368 369 END DO 369 370 END DO … … 393 394 & - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10)) 394 395 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) & 395 & * MAX(0.0,SIGN(1.0, nbot0(ji)-layer0+epsi20))396 & * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 396 397 END DO 397 398 END DO … … 440 441 DO jk = 1, nlay_s 441 442 DO ji = kideb, kiut 442 zswitch = MAX ( 0.0 , SIGN ( 1.0, epsi20- ht_s_b(ji) ) )443 zswitch = MAX ( 0.0 , SIGN ( 1.0, - ht_s_b(ji) ) ) 443 444 t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 ) 444 445 END DO … … 479 480 limsum = ( (icsuswi(ji)*(icsuind(ji)+jk-1) + & 480 481 (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) 481 zm0(ji,jk)= icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) &482 + limsum* zh_i(ji)483 END DO 484 END DO 485 486 DO ji = kideb, kiut 487 zm0(ji,nbot0(ji)) = icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) + dh_i_bott(ji) &488 + zh_i(ji) * nlayi0489 zm0(ji,1) = snicswi(ji)*dh_snowice(ji) +(1-snicswi(ji))*zm0(ji,1)482 zm0(ji,jk)= REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) & 483 + REAL(limsum) * zh_i(ji) 484 END DO 485 END DO 486 487 DO ji = kideb, kiut 488 zm0(ji,nbot0(ji)) = REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) + dh_i_bott(ji) & 489 + zh_i(ji) * REAL(nlayi0) 490 zm0(ji,1) = REAL(snicswi(ji))*dh_snowice(ji) + REAL(1-snicswi(ji))*zm0(ji,1) 490 491 END DO 491 492 … … 520 521 !---------------------------- 521 522 DO ji = kideb, kiut 522 ztmelts = ( 1.0- icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice523 & + icboswi(ji)* (-tmut * s_i_new(ji) ) & ! case of forming ice523 ztmelts = REAL( 1 - icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice 524 & + REAL( icboswi(ji) ) * (-tmut * s_i_new(ji) ) & ! case of forming ice 524 525 & + rtt ! in Kelvin 525 526 … … 527 528 ztform = t_i_b(ji,nlay_i) 528 529 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) ztform = t_bo_b(ji) 529 qm0(ji,nbot0(ji)) = ( 1.0- icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice530 & + icboswi(ji) * rhoic * ( cpic*(ztmelts-ztform) & ! case of forming ice530 qm0(ji,nbot0(ji)) = REAL( 1 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice 531 & + REAL( icboswi(ji) ) * rhoic * ( cpic*(ztmelts-ztform) & ! case of forming ice 531 532 + lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) ) & 532 533 - rcp*(ztmelts-rtt) ) * zthick0(ji,nbot0(ji) ) … … 539 540 ! energy of the flooding seawater 540 541 zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * & 541 (rhoic - rhosn) / rhoic * snicswi(ji) ! generally positive542 (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive 542 543 ! Heat conservation diagnostic 543 544 qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic … … 548 549 ! = enthalpy of snow + enthalpy of frozen water 549 550 zqsnic = zqsnow(ji) + zqsnic 550 qm0(ji,1) = snicswi(ji) * zqsnic + ( 1 - snicswi(ji) ) * qm0(ji,1) 551 qm0(ji,1) = REAL(snicswi(ji)) * zqsnic + REAL( 1 - snicswi(ji) ) * qm0(ji,1) 552 553 zji = MOD( npb(ji) - 1, jpi ) + 1 554 zjj = ( npb(ji) - 1 ) / jpi + 1 555 IF ( (zji.eq.jiindx).AND.(zjj.eq.jjindx) ) THEN 556 !clemclem 557 WRITE(numout,*) 'lim_thd_ent : qldif = ', qldif_1d(ji) 558 !clemclem 559 ENDIF 551 560 552 561 END DO ! ji … … 555 564 DO ji = kideb, kiut 556 565 ! Heat conservation 557 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06 +epsi20) ) &558 & * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + epsi20) )566 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06) ) & 567 & * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) 559 568 END DO 560 569 END DO … … 574 583 !------------------ 575 584 DO ji = kideb, kiut 576 zh_i(ji) = ht_i_b(ji) / nlay_i585 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 577 586 ENDDO 578 587 … … 605 614 q_i_b(ji,layer1) = q_i_b(ji,layer1) & 606 615 + zrl01(layer1,layer0)*qm0(ji,layer0) & 607 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06 +epsi20)) &608 * MAX(0.0,SIGN(1.0, nbot0(ji)-layer0+epsi20))616 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06)) & 617 * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 609 618 END DO 610 619 END DO -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r3294 r3938 29 29 USE lib_mpp ! MPP library 30 30 USE wrk_nemo ! work arrays 31 USE lib_fortran ! to use key_nosignedzero 31 32 32 33 IMPLICIT NONE … … 159 160 !Energy of melting q(S,T) [J.m-3] 160 161 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / & 161 MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * & 162 nlay_i 162 MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * REAL( nlay_i ) 163 163 zindb = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl))) !0 if no ice and 1 if yes 164 164 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl)*unit_fac*zindb … … 185 185 186 186 ! Default new ice thickness 187 DO jj = 1, jpj 188 DO ji = 1, jpi 189 hicol(ji,jj) = hiccrit(1) 190 END DO 191 END DO 187 hicol(:,:) = hiccrit(1) 192 188 193 189 IF (fraz_swi.eq.1.0) THEN … … 335 331 CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , & 336 332 jpi, jpj, npac(1:nbpac) ) 333 CALL tab_2d_1d( nbpac, rdmicif_1d (1:nbpac) , rdmicif, jpi, jpj, npac(1:nbpac) ) !martin 337 334 338 335 !------------------------------------------------------------------------------! … … 410 407 zdh_frazb(ji) = zfrazb*zv_newice(ji) 411 408 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 409 ! 410 !clem@mv rdmicif_1d(ji) = rdmicif_1d(ji) + zv_newice(ji) * rhoic !martin 412 411 END DO 413 412 … … 415 414 ! Salt flux due to new ice growth 416 415 !--------------------------------- 417 IF ( ( num_sal .EQ. 4 ) ) THEN 418 DO ji = 1, nbpac 419 zji = MOD( npac(ji) - 1, jpi ) + 1 420 zjj = ( npac(ji) - 1 ) / jpi + 1 421 fseqv_1d(ji) = fseqv_1d(ji) + & 422 ( sss_m(zji,zjj) - bulk_sal ) * rhoic * & 423 zv_newice(ji) / rdt_ice 424 END DO 425 ELSE 426 DO ji = 1, nbpac 427 zji = MOD( npac(ji) - 1, jpi ) + 1 428 zjj = ( npac(ji) - 1 ) / jpi + 1 429 fseqv_1d(ji) = fseqv_1d(ji) + & 430 ( sss_m(zji,zjj) - zs_newice(ji) ) * rhoic * & 431 zv_newice(ji) / rdt_ice 432 END DO ! ji 433 ENDIF 416 !clem@mv IF ( ( num_sal .EQ. 4 ) ) THEN 417 !clem@mv DO ji = 1, nbpac 418 !clem@mv ! IOVINO 419 !clem@mv fseqv_1d(ji) = fseqv_1d(ji) - bulk_sal * rhoic * & 420 !clem@mv zv_newice(ji) / rdt_ice 421 !clem@mv END DO 422 !clem@mv ELSE 423 !clem@mv DO ji = 1, nbpac 424 !clem@mv ! IOVINO 425 !clem@mv fseqv_1d(ji) = fseqv_1d(ji) - zs_newice(ji) * rhoic * & 426 !clem@mv zv_newice(ji) / rdt_ice 427 !clem@mv END DO ! ji 428 !clem@mv ENDIF 434 429 435 430 !------------------------------------ … … 460 455 zjj = ( npac(ji) - 1 ) / jpi + 1 461 456 diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice 462 END DO !ji457 END DO !ji 463 458 464 459 !------------------------------------------------------------------------------! … … 480 475 DO ji = 1, nbpac 481 476 ! vectorize 482 IF ( za_newice(ji) .GT. ( 1.0- zat_i_ac(ji) ) ) THEN483 zda_res(ji) = za_newice(ji) - ( 1.0- zat_i_ac(ji) )477 IF ( za_newice(ji) .GT. ( amax - zat_i_ac(ji) ) ) THEN 478 zda_res(ji) = za_newice(ji) - ( amax - zat_i_ac(ji) ) 484 479 zdv_res(ji) = zda_res(ji) * zh_newice(ji) 485 480 za_newice(ji) = za_newice(ji) - zda_res(ji) 486 481 zv_newice(ji) = zv_newice(ji) - zdv_res(ji) 487 482 ELSE 488 483 zda_res(ji) = 0.0 489 484 zdv_res(ji) = 0.0 … … 512 507 DO ji = 1, nbpac 513 508 jl = zcatac(ji) ! categroy in which new ice is put 514 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) ) ) ! zindb=1 if ice =0 otherwise509 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 515 510 zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness 516 511 zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness 517 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi1 1) ) ! ice totally new in jl category512 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! ice totally new in jl category 518 513 END DO 519 514 … … 522 517 jl = zcatac(ji) 523 518 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 524 zalphai = MIN( zhice_old(ji,jl) * jk / nlay_i, zh_newice(ji) ) &525 & - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i, zh_newice(ji) )519 zalphai = MIN( zhice_old(ji,jl) * REAL( jk ) / REAL( nlay_i ), zh_newice(ji) ) & 520 & - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) ) 526 521 ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) & 527 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / nlay_i&522 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / REAL( nlay_i ) & 528 523 + za_newice(ji) * ze_newice(ji) * zalphai & 529 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( ( zv_i_ac(ji,jl) ) / nlay_i)524 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) ) 530 525 END DO 531 526 END DO … … 563 558 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 564 559 DO ji = 1, nbpac 565 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) ) ) ! zindb=1 if ice =0 otherwise560 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 566 561 zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 567 562 zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb & … … 575 570 DO jk = 1, nlay_i 576 571 DO ji = 1, nbpac 577 zthick0(ji,jk,jl) = zhice_old(ji,jl) / nlay_i572 zthick0(ji,jk,jl) = zhice_old(ji,jl) / REAL( nlay_i ) 578 573 zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl) 579 574 END DO … … 594 589 DO layer = 1, nlay_i + 1 595 590 DO ji = 1, nbpac 596 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) )591 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 597 592 ! Redistributing energy on the new grid 598 zweight = MAX ( MIN( zhice_old(ji,jl) * layer , zdummy(ji,jl) * jk) &599 & - MAX( zhice_old(ji,jl) * ( layer - 1 ) , zdummy(ji,jl) *( jk - 1 ) ) , 0._wp ) &600 & /( MAX( nlay_i* zthick0(ji,layer,jl),epsi10) ) * zindb593 zweight = MAX ( MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) ) & 594 & - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp ) & 595 & /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb 601 596 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 602 597 END DO ! ji … … 608 603 DO jk = 1, nlay_i 609 604 DO ji = 1, nbpac 610 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )605 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 611 606 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) & 612 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * nlay_i* zindb607 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb 613 608 END DO 614 609 END DO … … 620 615 DO jl = 1, jpl 621 616 DO ji = 1, nbpac 622 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes617 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 623 618 zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 624 619 END DO … … 631 626 DO jl = 1, jpl 632 627 DO ji = 1, nbpac 633 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes628 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 634 629 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 635 630 zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb … … 637 632 END DO 638 633 ENDIF 634 635 !-------------------------------- 636 ! Update mass/salt fluxes (clem) 637 !-------------------------------- 638 DO jl = 1, jpl 639 DO ji = 1, nbpac 640 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 641 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 642 rdmicif_1d(ji) = rdmicif_1d(ji) + zdv * rhoic !* zindb 643 fseqv_1d(ji) = fseqv_1d(ji) - zdv * rhoic * zs_newice(ji) / rdt_ice * zindb 644 END DO 645 END DO 639 646 640 647 !------------------------------------------------------------------------------! … … 652 659 END DO 653 660 CALL tab_1d_2d( nbpac, fseqv , npac(1:nbpac), fseqv_1d (1:nbpac) , jpi, jpj ) 661 CALL tab_1d_2d( nbpac, rdmicif , npac(1:nbpac), rdmicif_1d (1:nbpac) , jpi, jpj ) 654 662 ! 655 663 ENDIF ! nbpac > 0 … … 660 668 DO jl = 1, jpl 661 669 DO jk = 1, nlay_i ! heat content in 10^9 Joules 662 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i/ unit_fac670 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac 663 671 END DO 664 672 END DO -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
r3294 r3938 24 24 USE lib_mpp ! MPP library 25 25 USE wrk_nemo ! work arrays 26 USE lib_fortran ! to use key_nosignedzero 26 27 27 28 IMPLICIT NONE … … 31 32 PUBLIC lim_thd_sal_init ! called by iceini module 32 33 34 REAL(wp) :: epsi20 = 1e-20_wp ! constant values 33 35 !!---------------------------------------------------------------------- 34 36 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 67 69 IF( num_sal == 1 ) THEN 68 70 ! 69 DO jk = 1, nlay_i 70 DO ji = kideb, kiut 71 s_i_b(ji,jk) = bulk_sal 72 END DO ! ji 73 END DO ! jk 74 ! 75 DO ji = kideb, kiut 76 sm_i_b(ji) = bulk_sal 77 END DO ! ji 78 ! 71 s_i_b (kideb:kiut,1:nlay_i) = bulk_sal 72 sm_i_b (kideb:kiut) = bulk_sal 73 s_i_new(kideb:kiut) = bulk_sal 74 ! 79 75 ENDIF 80 76 … … 90 86 DO ji = kideb, kiut 91 87 zhiold(ji) = ht_i_b(ji) - dh_i_bott(ji) - dh_snowice(ji) - dh_i_surf(ji) 88 zsiold(ji) = sm_i_b(ji) 92 89 END DO 93 90 … … 98 95 DO jk = 1, nlay_i 99 96 DO ji = kideb, kiut 100 ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i97 ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / REAL (nlay_i ) 101 98 END DO 102 99 END DO … … 125 122 ! only drainage terms ( gravity drainage and flushing ) 126 123 ! snow ice / bottom sources are added in lim_thd_ent to conserve energy 127 zsiold(ji) = sm_i_b(ji)128 124 sm_i_b(ji) = sm_i_b(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji) 129 125 … … 156 152 ! Salt flux - brine drainage 157 153 !---------------------------- 158 DO ji = kideb, kiut154 DO ji = kideb, kiut 159 155 i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) ) 160 fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) & 161 & * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice 156 fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice 157 !i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - zhiold(ji) ) ) 158 !fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * zhiold(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice 159 !clem fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) & 160 !clem & * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice 162 161 IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp 163 END DO ! ji162 END DO ! ji 164 163 165 164 ! Only necessary for conservation check since salinity is modified … … 211 210 ENDIF ! num_sal 212 211 213 !------------------------------------------------------------------------------|214 ! 5) Computation of salt flux due to Bottom growth215 !------------------------------------------------------------------------------|216 217 IF ( num_sal == 4 ) THEN218 DO ji = kideb, kiut219 zji = MOD( npb(ji) - 1 , jpi ) + 1220 zjj = ( npb(ji) - 1 ) / jpi + 1221 fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal ) &222 & * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice223 END DO224 ELSE225 DO ji = kideb, kiut226 zji = MOD( npb(ji) - 1 , jpi ) + 1227 zjj = ( npb(ji) - 1 ) / jpi + 1228 fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - s_i_new(ji) ) &229 & * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice230 END DO231 ENDIF232 212 ! 233 213 CALL wrk_dealloc( jpij, ze_init, zhiold, zsiold ) -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r3294 r3938 27 27 USE wrk_nemo ! work arrays 28 28 USE prtctl ! Print control 29 USE lib_fortran ! to use key_nosignedzero 30 USE limvar ! clem for ice thickness correction 29 31 30 32 IMPLICIT NONE … … 35 37 REAL(wp) :: epsi06 = 1.e-06_wp ! constant values 36 38 REAL(wp) :: epsi03 = 1.e-03_wp 37 REAL(wp) :: zeps10 = 1.e-10_wp39 REAL(wp) :: epsi10 = 1.e-10_wp 38 40 REAL(wp) :: epsi16 = 1.e-16_wp 41 REAL(wp) :: epsi20 = 1.e-20_wp 39 42 REAL(wp) :: rzero = 0._wp 40 43 REAL(wp) :: rone = 1._wp … … 65 68 INTEGER, INTENT(in) :: kt ! number of iteration 66 69 ! 67 INTEGER :: ji, jj, jk, jl, layer ! dummy loop indices 70 INTEGER :: ji, jj, jk, jl, jm, layer ! dummy loop indices 71 INTEGER :: jbnd1, jbnd2 68 72 INTEGER :: initad ! number of sub-timestep for the advection 69 73 INTEGER :: ierr ! error status 70 REAL(wp) :: zindb , zindsn , zindic ! local scalar74 REAL(wp) :: zindb , zindsn , zindic, zindh, zinda ! local scalar 71 75 REAL(wp) :: zusvosn, zusvoic, zbigval ! - - 72 76 REAL(wp) :: zcfl , zusnit , zrtt ! - - … … 76 80 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 77 81 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e 82 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 83 ! mass and salt flux (clem) 84 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold ! old ice volume... 85 ! correct ice thickness (clem) 86 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration 87 REAL(wp) :: zdv, zda, zvi, zvs, zsmv 78 88 !!--------------------------------------------------------------------- 79 89 … … 81 91 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 82 92 CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e ) 93 94 CALL wrk_alloc( jpi,jpj,jpl,zviold ) ! clem 95 CALL wrk_alloc( jpi,jpj,jpl,zaiold, zhimax ) ! clem 96 97 ! ------------------------------- 98 !- check conservation (C Rousset) 99 IF (ln_limdiahsb) THEN 100 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 101 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 102 zchk_fw_b = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) 103 zchk_fs_b = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) 104 ENDIF 105 !- check conservation (C Rousset) 106 ! ------------------------------- 83 107 84 108 IF( numit == nstart .AND. lwp ) THEN … … 95 119 IF( ln_limdyn ) THEN ! Advection of sea ice properties ! 96 120 ! !-------------------------------------! 97 ! 98 121 ! mass and salt flux init (clem) 122 zviold(:,:,:) = v_i(:,:,:) 123 124 !--- Thickness correction init. (clem) ------------------------------- 125 CALL lim_var_glo2eqv 126 zaiold(:,:,:) = a_i(:,:,:) 127 !--------------------------------------------------------------------- 128 ! Record max of the surrounding ice thicknesses for correction in limupdate 129 ! in case advection creates ice too thick. 130 !--------------------------------------------------------------------- 131 zhimax(:,:,:) = ht_i(:,:,:) 132 DO jl = 1, jpl 133 DO jj = 2, jpjm1 134 DO ji = 2, jpim1 135 zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) ) 136 END DO 137 END DO 138 CALL lbc_lnk(zhimax(:,:,jl),'T',1.) 139 END DO 140 99 141 !------------------------- 100 142 ! transported fields … … 125 167 ! ENDIF 126 168 !!gm end 127 initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )169 initad = 1 + NINT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) ) 128 170 zusnit = 1.0 / REAL( initad ) 129 171 IF( zcfl > 0.5 .AND. lwp ) & … … 134 176 DO jk = 1,initad 135 177 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 136 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )178 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 137 179 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), & 138 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )180 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 139 181 DO jl = 1, jpl 140 182 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- … … 174 216 ELSE 175 217 DO jk = 1, initad 176 CALL lim_adv_y( zusnit, v_ice, r zero, zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area177 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )178 CALL lim_adv_x( zusnit, u_ice, r one, zsm, zs0ow (:,:), sxopw(:,:), &179 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )218 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 219 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 220 CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0ow (:,:), sxopw(:,:), & 221 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 180 222 DO jl = 1, jpl 181 CALL lim_adv_y( zusnit, v_ice, r zero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---223 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 182 224 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 183 CALL lim_adv_x( zusnit, u_ice, r one, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &225 CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 184 226 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 185 CALL lim_adv_y( zusnit, v_ice, r zero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---227 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 186 228 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 187 CALL lim_adv_x( zusnit, u_ice, r one, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &229 CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 188 230 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 189 CALL lim_adv_y( zusnit, v_ice, r zero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---231 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 190 232 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 191 CALL lim_adv_x( zusnit, u_ice, r one, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &233 CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 192 234 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 193 194 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 235 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 195 236 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 196 CALL lim_adv_x( zusnit, u_ice, r one, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &237 CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 197 238 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 198 CALL lim_adv_y( zusnit, v_ice, r zero, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---239 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 199 240 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 200 CALL lim_adv_x( zusnit, u_ice, r one, zsm, zs0a (:,:,jl), sxa (:,:,jl), &241 CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0a (:,:,jl), sxa (:,:,jl), & 201 242 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 202 CALL lim_adv_y( zusnit, v_ice, r zero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---243 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 203 244 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 204 CALL lim_adv_x( zusnit, u_ice, r one, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &245 CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 205 246 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 206 247 DO layer = 1, nlay_i !--- ice heat contents --- 207 CALL lim_adv_y( zusnit, v_ice, r zero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &248 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 208 249 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 209 250 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 210 CALL lim_adv_x( zusnit, u_ice, r one, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &251 CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 211 252 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 212 253 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) … … 270 311 END DO 271 312 272 CALL lim_hdf( zs0ice (:,:,jl) )273 CALL lim_hdf( zs0sn (:,:,jl) )274 CALL lim_hdf( zs0sm (:,:,jl) )275 CALL lim_hdf( zs0oi (:,:,jl) )276 CALL lim_hdf( zs0a (:,:,jl) )277 CALL lim_hdf( zs0c0 (:,:,jl) )278 DO jk = 1, nlay_i279 CALL lim_hdf( zs0e (:,:,jk,jl) )280 END DO313 CALL lim_hdf( zs0ice (:,:,jl) ) 314 CALL lim_hdf( zs0sn (:,:,jl) ) 315 CALL lim_hdf( zs0sm (:,:,jl) ) 316 CALL lim_hdf( zs0oi (:,:,jl) ) 317 CALL lim_hdf( zs0a (:,:,jl) ) 318 CALL lim_hdf( zs0c0 (:,:,jl) ) 319 DO jk = 1, nlay_i 320 CALL lim_hdf( zs0e (:,:,jk,jl) ) 321 END DO 281 322 END DO 282 323 … … 284 325 ! Remultiply everything by ice area 285 326 !----------------------------------------- 286 zs0ow(:,:) = MAX( rzero, zs0ow(:,:) * area(:,:) )287 DO jl = 1, jpl288 zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) ) !!bug: est-ce utile289 zs0sn (:,:,jl) = MAX( rzero, zs0sn (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres290 zs0sm (:,:,jl) = MAX( rzero, zs0sm (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres291 zs0oi (:,:,jl) = MAX( rzero, zs0oi (:,:,jl) * area(:,:) )292 zs0a (:,:,jl) = MAX( rzero, zs0a (:,:,jl) * area(:,:) ) !! suppress both change le resultat293 zs0c0 (:,:,jl) = MAX( rzero, zs0c0 (:,:,jl) * area(:,:) )294 DO jk = 1, nlay_i295 zs0e(:,:,jk,jl) = MAX( rzero, zs0e (:,:,jk,jl) * area(:,:) )296 END DO ! jk297 END DO ! jl327 !clem zs0ow(:,:) = MAX( rzero, zs0ow(:,:) * area(:,:) ) 328 !clem DO jl = 1, jpl 329 !clem zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) ) !!bug: est-ce utile 330 !clem zs0sn (:,:,jl) = MAX( rzero, zs0sn (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres 331 !clem zs0sm (:,:,jl) = MAX( rzero, zs0sm (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres 332 !clem zs0oi (:,:,jl) = MAX( rzero, zs0oi (:,:,jl) * area(:,:) ) 333 !clem zs0a (:,:,jl) = MAX( rzero, zs0a (:,:,jl) * area(:,:) ) !! suppress both change le resultat 334 !clem zs0c0 (:,:,jl) = MAX( rzero, zs0c0 (:,:,jl) * area(:,:) ) 335 !clem DO jk = 1, nlay_i 336 !clem zs0e(:,:,jk,jl) = MAX( rzero, zs0e (:,:,jk,jl) * area(:,:) ) 337 !clem END DO ! jk 338 !clem END DO ! jl 298 339 299 340 !------------------------------------------------------------------------------! … … 305 346 !-------------------------------------------------- 306 347 307 DO jl = 1, jpl308 DO jk = 1, nlay_i309 DO jj = 1, jpj310 DO ji = 1, jpi311 zs0e(ji,jj,jk,jl) = MAX( rzero, zs0e(ji,jj,jk,jl) / area(ji,jj) )312 END DO313 END DO314 END DO315 END DO316 317 DO jj = 1, jpj318 DO ji = 1, jpi319 zs0ow(ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) )320 END DO321 END DO348 !clem DO jl = 1, jpl 349 !clem DO jk = 1, nlay_i 350 !clem DO jj = 1, jpj 351 !clem DO ji = 1, jpi 352 !clem zs0e(ji,jj,jk,jl) = MAX( rzero, zs0e(ji,jj,jk,jl) / area(ji,jj) ) 353 !clem END DO 354 !clem END DO 355 !clem END DO 356 !clem END DO 357 358 !clem DO jj = 1, jpj 359 !clem DO ji = 1, jpi 360 !clem zs0ow(ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) ) 361 !clem END DO 362 !clem END DO 322 363 323 364 zs0at(:,:) = 0._wp … … 325 366 DO jj = 1, jpj 326 367 DO ji = 1, jpi 327 zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl)/area(ji,jj) ) 328 zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl)/area(ji,jj) ) 329 zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl)/area(ji,jj) ) 330 zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl)/area(ji,jj) ) 331 zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl)/area(ji,jj) ) 332 zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl)/area(ji,jj) ) 368 !clem zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl)/area(ji,jj) ) 369 !clem zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl)/area(ji,jj) ) 370 !clem zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl)/area(ji,jj) ) 371 !clem zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl)/area(ji,jj) ) 372 !clem zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl)/area(ji,jj) ) 373 !clem zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl)/area(ji,jj) ) 374 zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl) ) 375 zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl) ) 376 zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl) ) 377 zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl) ) 378 zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl) ) 379 zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl) ) 333 380 zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) 334 381 END DO … … 341 388 DO jj = 1, jpj 342 389 DO ji = 1, jpi 343 zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - zeps10) )390 zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - epsi10) ) 344 391 zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp ) 345 392 ato_i(ji,jj) = zs0ow(ji,jj) … … 347 394 END DO 348 395 396 ! 397 ! 349 398 DO jl = 1, jpl ! Remove very small areas 350 399 DO jj = 1, jpj 351 400 DO ji = 1, jpi 352 zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - zeps10) ) 401 zvi = zs0ice(ji,jj,jl) 402 zvs = zs0sn(ji,jj,jl) 353 403 ! 354 zs0a(ji,jj,jl) = zindb * MIN( zs0a(ji,jj,jl), 0.99 ) 404 zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - epsi10) ) 405 ! 406 !zs0a(ji,jj,jl) = zindb * MIN( zs0a(ji,jj,jl), 0.99 ) 355 407 v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl) 356 408 v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl) 357 409 ! 358 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )359 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )410 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 411 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 360 412 zindb = MAX( zindsn, zindic ) 413 ! 361 414 zs0a(ji,jj,jl) = zindb * zs0a(ji,jj,jl) !ice concentration 362 415 a_i (ji,jj,jl) = zs0a(ji,jj,jl) 363 416 v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl) 364 417 v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl) 418 ! 419 ! Update mass fluxes (clem) 420 rdmicif(ji,jj) = rdmicif(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 421 rdmsnif(ji,jj) = rdmsnif(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 422 END DO 423 END DO 424 END DO 425 426 !--- Thickness correction in case too high (clem) -------------------------------------------------------- 427 CALL lim_var_glo2eqv 428 DO jl = 1, jpl 429 DO jj = 1, jpj 430 DO ji = 1, jpi 431 432 IF ( v_i(ji,jj,jl) > 0. ) THEN 433 zvi = v_i(ji,jj,jl) 434 zvs = v_s(ji,jj,jl) 435 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 436 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 437 438 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 439 ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 440 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 441 zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 442 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 ) 443 ELSE 444 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 445 zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 446 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 ) 447 ENDIF 448 449 !zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 450 v_i(ji,jj,jl) = a_i(ji,jj,jl) * ht_i(ji,jj,jl) 451 v_s(ji,jj,jl) = a_i(ji,jj,jl) * ht_s(ji,jj,jl) 452 453 ! Update mass fluxes (clem) 454 rdmicif(ji,jj) = rdmicif(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 455 rdmsnif(ji,jj) = rdmsnif(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 456 457 ENDIF 458 459 diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) / rdt_ice 460 365 461 END DO 366 462 END DO 367 463 END DO 368 464 465 ! --- 369 466 DO jj = 1, jpj 370 467 DO ji = 1, jpi 371 zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) 468 ! ato_i(ji,jj) = 1._wp - SUM( a_i(ji,jj,1:jpl) ) !clem@rm-ow-advection 469 zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) ! clem@useless?? 372 470 END DO 373 471 END DO … … 377 475 !---------------------- 378 476 379 zbigval = 1. d+13477 zbigval = 1.e+13 380 478 381 479 DO jl = 1, jpl 382 480 DO jj = 1, jpj 383 481 DO ji = 1, jpi 482 zsmv = zs0sm(ji,jj,jl) 384 483 385 484 ! Switches and dummy variables … … 387 486 zusvoic = 1.0/MAX( v_i(ji,jj,jl) , epsi16 ) 388 487 zrtt = 173.15 * rone 389 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )390 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )488 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 489 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 391 490 zindb = MAX( zindsn, zindic ) 392 491 … … 394 493 zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj) , & 395 494 zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl) 396 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 397 smv_i(ji,jj,jl) = zindic*zsal + (1.0-zindic)*0.0 495 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) smv_i(ji,jj,jl) = zindic*zsal 398 496 399 497 zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / & … … 403 501 ! Snow heat content 404 502 ze = MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval ) 405 e_s(ji,jj,1,jl) = zindsn * ze + (1.0 - zindsn) * 0.0 406 503 e_s(ji,jj,1,jl) = zindsn * ze 504 505 ! Update salt fluxes (clem) 506 fsalt_res(ji,jj) = fsalt_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic / rdt_ice 407 507 END DO !ji 408 508 END DO !jj … … 414 514 DO ji = 1, jpi 415 515 ! Ice heat content 416 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )516 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 417 517 ze = MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval ) 418 e_i(ji,jj,jk,jl) = zindic * ze + ( 1.0 - zindic ) * 0.0518 e_i(ji,jj,jk,jl) = zindic * ze 419 519 END DO !ji 420 520 END DO ! jj 421 521 END DO ! jk 422 522 END DO ! jl 523 524 525 ! --- agglomerate variables (clem) ----------------- 526 vt_i (:,:) = 0._wp 527 vt_s (:,:) = 0._wp 528 at_i (:,:) = 0._wp 529 ! 530 DO jl = 1, jpl 531 DO jj = 1, jpj 532 DO ji = 1, jpi 533 ! 534 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 535 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 536 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 537 ! 538 zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi16 ) ) 539 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda ! ice thickness 540 END DO 541 END DO 542 END DO 543 ! ------------------------------------------------- 544 545 423 546 424 547 ENDIF … … 455 578 END DO 456 579 ENDIF 580 ! ------------------------------- 581 !- check conservation (C Rousset) 582 IF (ln_limdiahsb) THEN 583 zchk_fs = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 584 zchk_fw = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 585 586 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 587 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 588 589 IF(lwp) THEN 590 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limtrp) = ',(zchk_v_i * 86400.) 591 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * 86400.) 592 IF ( MINVAL( v_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limtrp) = ',(MINVAL(v_i) * 1.e-3) 593 IF ( MINVAL( a_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation a_i<0 (limtrp) = ',MINVAL(a_i) 594 ENDIF 595 ENDIF 596 !- check conservation (C Rousset) 597 ! ------------------------------- 457 598 ! 458 599 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow ) 459 600 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 460 601 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e ) 602 603 CALL wrk_dealloc( jpi,jpj,jpl,zaiold, zhimax ) ! clem 461 604 ! 462 605 END SUBROUTINE lim_trp -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r3294 r3938 53 53 USE lib_mpp ! MPP library 54 54 USE wrk_nemo ! work arrays 55 USE lib_fortran ! to use key_nosignedzero 55 56 56 57 IMPLICIT NONE … … 61 62 PUBLIC lim_var_eqv2glo ! 62 63 PUBLIC lim_var_salprof ! 64 PUBLIC lim_var_icetm ! 63 65 PUBLIC lim_var_bv ! 64 66 PUBLIC lim_var_salprof1d ! 65 67 66 REAL(wp) :: eps 20 = 1.e-20_wp ! module constants67 REAL(wp) :: eps 16 = 1.e-16_wp ! - -68 REAL(wp) :: eps 13 = 1.e-13_wp ! - -69 REAL(wp) :: eps 10 = 1.e-10_wp ! - -70 REAL(wp) :: eps 06 = 1.e-06_wp ! - -68 REAL(wp) :: epsi20 = 1.e-20_wp ! module constants 69 REAL(wp) :: epsi16 = 1.e-16_wp ! - - 70 REAL(wp) :: epsi13 = 1.e-13_wp ! - - 71 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 72 REAL(wp) :: epsi06 = 1.e-06_wp ! - - 71 73 REAL(wp) :: zzero = 0.e0 ! - - 72 74 REAL(wp) :: zone = 1.e0 ! - - … … 96 98 ! 97 99 INTEGER :: ji, jj, jk, jl ! dummy loop indices 98 REAL(wp) :: zinda 100 REAL(wp) :: zinda, zindb 99 101 !!------------------------------------------------------------------ 100 102 … … 115 117 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 116 118 ! 117 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10) )118 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , eps 16 ) * zinda ! ice thickness119 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) ) 120 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda ! ice thickness 119 121 END DO 120 122 END DO … … 136 138 DO jj = 1, jpj 137 139 DO ji = 1, jpi 140 zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi16 ) ) 141 zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) ) 138 142 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 139 zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - 0.10 ) ) 140 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , eps13 ) * zinda ! ice salinity 141 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) ) 142 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , eps13 ) * zinda ! ice age 143 END DO 144 END DO 145 END DO 146 ! 143 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi16 ) * zinda ! ice salinity 144 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi16 ) * zindb ! ice age 145 END DO 146 END DO 147 END DO 148 ! 147 149 DO jl = 1, jpl 148 150 DO jk = 1, nlay_i … … 174 176 DO jj = 1, jpj 175 177 DO ji = 1, jpi 176 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) ) ) !0 if no ice and 1 if yes177 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps 10 ) * zindb178 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps 10 ) * zindb179 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps 10 ) * zindb178 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 179 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb 180 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb 181 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb 180 182 END DO 181 183 END DO … … 186 188 DO jj = 1, jpj 187 189 DO ji = 1, jpi 188 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) ) ) !0 if no ice and 1 if yes189 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , eps 10 ) * zindb190 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 191 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * zindb 190 192 END DO 191 193 END DO … … 207 209 DO ji = 1, jpi 208 210 ! ! Energy of melting q(S,T) [J.m-3] 209 zq_i = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , eps 06 ) * REAL(nlay_i,wp)211 zq_i = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi06 ) * REAL(nlay_i,wp) 210 212 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) ) ! zindb = 0 if no ice and 1 if yes 211 213 zq_i = zq_i * unit_fac * zindb !convert units … … 233 235 DO ji = 1, jpi 234 236 !Energy of melting q(S,T) [J.m-3] 235 zq_s = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , eps 06 ) ) * REAL(nlay_s,wp)237 zq_s = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL(nlay_s,wp) 236 238 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) ) ) ! zindb = 0 if no ice and 1 if yes 237 239 zq_s = zq_s * unit_fac * zindb ! convert units … … 252 254 DO jj = 1, jpj 253 255 DO ji = 1, jpi 254 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , -a_i(ji,jj,jl) ) ) ) & 255 & * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) ) ) ) 256 tm_i(ji,jj) = tm_i(ji,jj) + t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 257 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , eps10 ) ) 256 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 257 tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 258 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 258 259 END DO 259 260 END DO … … 337 338 DO ji = 1, jpi 338 339 ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise 339 zind0 = MAX( 0. 0 , SIGN( 1.0, s_i_0 - sm_i(ji,jj,jl) ) )340 zind0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i(ji,jj,jl) ) ) 340 341 ! zind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws 341 zind01 = ( 1. 0 - zind0 ) * MAX( 0.0 , SIGN( 1.0, s_i_1 - sm_i(ji,jj,jl) ) )342 zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i(ji,jj,jl) ) ) 342 343 ! If 2.sm_i GE sss_m then zindbal = 1 343 zindbal = MAX( 0. 0 , SIGN( 1.0 , 2.* sm_i(ji,jj,jl) - sss_m(ji,jj) ) )344 zalpha(ji,jj,jl) = zind0 * 1.0+ zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )345 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1. 0- zindbal )346 END DO 347 END DO 348 END DO 349 350 dummy_fac = 1._wp / nlay_i! Computation of the profile344 zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 345 zalpha(ji,jj,jl) = zind0 + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 ) 346 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zindbal ) 347 END DO 348 END DO 349 END DO 350 351 dummy_fac = 1._wp / REAL( nlay_i ) ! Computation of the profile 351 352 DO jl = 1, jpl 352 353 DO jk = 1, nlay_i … … 388 389 389 390 391 SUBROUTINE lim_var_icetm 392 !!------------------------------------------------------------------ 393 !! *** ROUTINE lim_var_icetm *** 394 !! 395 !! ** Purpose : computes mean sea ice temperature 396 !!------------------------------------------------------------------ 397 INTEGER :: ji, jj, jk, jl ! dummy loop indices 398 REAL(wp) :: zindb ! - - 399 !!------------------------------------------------------------------ 400 401 ! Mean sea ice temperature 402 tm_i(:,:) = 0._wp 403 DO jl = 1, jpl 404 DO jk = 1, nlay_i 405 DO jj = 1, jpj 406 DO ji = 1, jpi 407 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 408 tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 409 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 410 END DO 411 END DO 412 END DO 413 END DO 414 415 END SUBROUTINE lim_var_icetm 416 417 390 418 SUBROUTINE lim_var_bv 391 419 !!------------------------------------------------------------------ … … 399 427 !!------------------------------------------------------------------ 400 428 INTEGER :: ji, jj, jk, jl ! dummy loop indices 401 REAL(wp) :: zbvi, zind b ! local scalars429 REAL(wp) :: zbvi, zinda, zindb ! local scalars 402 430 !!------------------------------------------------------------------ 403 431 ! … … 407 435 DO jj = 1, jpj 408 436 DO ji = 1, jpi 409 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl))) !0 if no ice and 1 if yes 410 zbvi = - zindb * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - 273.15 , eps13 ) & 437 zinda = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi16 ) ) ) 438 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi16 ) ) ) 439 zbvi = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi16 ) & 411 440 & * v_i(ji,jj,jl) / REAL(nlay_i,wp) 412 bv_i(ji,jj) = bv_i(ji,jj) + z bvi / MAX( vt_i(ji,jj) , eps13)441 bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi / MAX( vt_i(ji,jj) , epsi16 ) 413 442 END DO 414 443 END DO -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r3294 r3938 10 10 !! lim_wri : write of the diagnostics variables in ouput file 11 11 !! lim_wri_init : initialization and namelist read 12 !! lim_wri_state : write for initial state or/and abandon 12 13 !!---------------------------------------------------------------------- 13 14 USE ioipsl … … 25 26 USE wrk_nemo ! work arrays 26 27 USE par_ice 28 USE iom 29 USE timing ! Timing 30 USE lib_fortran ! Fortran utilities 27 31 28 32 IMPLICIT NONE … … 30 34 31 35 PUBLIC lim_wri ! routine called by lim_step.F90 32 33 INTEGER, PARAMETER :: jpnoumax = 40 !: maximum number of variable for ice output 36 PUBLIC lim_wri_state ! called by dia_wri_state 37 38 INTEGER, PARAMETER :: jpnoumax = 43 !: maximum number of variable for ice output 34 39 35 40 INTEGER :: noumef ! number of fields … … 47 52 INTEGER , DIMENSION(jpnoumax) :: nc , nca ! switch for saving field ( = 1 ) or not ( = 0 ) 48 53 49 REAL(wp) :: epsi 16 = 1e-16_wp54 REAL(wp) :: epsi06 = 1e-6_wp 50 55 REAL(wp) :: zzero = 0._wp 51 56 REAL(wp) :: zone = 1._wp … … 76 81 INTEGER :: ierr 77 82 REAL(wp),DIMENSION(1) :: zdept 78 REAL(wp) :: zsto, zjulian, zout, zindh, zinda, zindb 83 REAL(wp) :: zsto, zjulian, zout, zindh, zinda, zindb, zindc 79 84 REAL(wp), POINTER, DIMENSION(:,:,:) :: zcmo, zcmoa 80 85 REAL(wp), POINTER, DIMENSION(:,: ) :: zfield … … 88 93 INTEGER , ALLOCATABLE, DIMENSION(:), SAVE :: ndexitd 89 94 !!------------------------------------------------------------------- 95 96 IF( nn_timing == 1 ) CALL timing_start('limwri') 90 97 91 98 CALL wrk_alloc( jpi, jpj, zfield ) … … 126 133 CALL ymds2ju ( nyear, nmonth, nday, rdt, zjulian ) 127 134 zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment 128 CALL dia_nam ( clhstnam, nwrite, 'icemod ' )135 CALL dia_nam ( clhstnam, nwrite, 'icemod_old' ) 129 136 CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit, 1, jpi, 1, jpj, niter, zjulian, rdt_ice, & 130 137 & nhorid, nice, domain_id=nidom, snc4chunks=snc4set ) … … 158 165 nhorida, & ! ? linked with horizontal ... 159 166 nicea , domain_id=nidom, snc4chunks=snc4set) ! file 160 CALL histvert( nicea, "icethi", "L levels", & 161 "m", ipl , hi_mean , nz ) 167 CALL histvert( nicea, "icethi", "L levels","m", ipl , hi_mean , nz ) 162 168 DO jl = 1, jpl 163 169 zmaskitd(:,:,jl) = tmask(:,:,1) … … 197 203 zcmoa( 1:jpi, 1:jpj, 1:jpnoumax ) = 0._wp 198 204 205 ! Ice surface temperature and some fluxes 199 206 DO jl = 1, jpl 200 207 DO jj = 1, jpj 201 208 DO ji = 1, jpi 202 zindh = MAX( zzero , SIGN( zone , vt_i(ji,jj) * at_i(ji,jj) - 0.10 ) ) 203 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) ) 209 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi06 ) ) 204 210 zcmo(ji,jj,17) = zcmo(ji,jj,17) + a_i(ji,jj,jl)*qsr_ice (ji,jj,jl) 205 211 zcmo(ji,jj,18) = zcmo(ji,jj,18) + a_i(ji,jj,jl)*qns_ice(ji,jj,jl) 206 zcmo(ji,jj,27) = zcmo(ji,jj,27) + t_su(ji,jj,jl)*a_i(ji,jj,jl)/MAX(at_i(ji,jj),epsi16)*zinda 212 zcmo(ji,jj,27) = zcmo(ji,jj,27) + zinda*(t_su(ji,jj,jl)-rtt)*a_i(ji,jj,jl)/MAX(at_i(ji,jj),epsi06) 213 zcmo(ji,jj,21) = zcmo(ji,jj,21) + zinda*oa_i(ji,jj,jl)/MAX(at_i(ji,jj),epsi06) 207 214 END DO 208 215 END DO 209 216 END DO 210 217 218 ! Mean sea ice temperature 219 CALL lim_var_icetm 220 221 ! Brine volume 211 222 CALL lim_var_bv 212 223 213 224 DO jj = 2 , jpjm1 214 225 DO ji = 2 , jpim1 215 zindh = MAX( zzero , SIGN( zone , vt_i(ji,jj) * at_i(ji,jj) - 0.10 ) ) 216 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) ) 217 zindb = zindh * zinda 226 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi06 ) ) 227 zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) ) ) 218 228 219 229 zcmo(ji,jj,1) = at_i(ji,jj) 220 zcmo(ji,jj,2) = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi 16 ) * zinda221 zcmo(ji,jj,3) = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi 16 ) * zinda222 zcmo(ji,jj,4) = diag_bot_gr(ji,jj) * 86400.0 * zinda! Bottom thermodynamic ice production223 zcmo(ji,jj,5) = diag_dyn_gr(ji,jj) * 86400.0 * zinda! Dynamic ice production (rid/raft)224 zcmo(ji,jj,22) = diag_lat_gr(ji,jj) * 86400.0 * zinda! Lateral thermodynamic ice production225 zcmo(ji,jj,23) = diag_sni_gr(ji,jj) * 86400.0 * zinda! Snow ice production ice production226 zcmo(ji,jj,24) = tm_i(ji,jj) - rtt227 228 zcmo(ji,jj,6) = fbif 229 zcmo(ji,jj,7) = zindb *( u_ice(ji,jj) * tmu(ji,jj) + u_ice(ji-1,jj) * tmu(ji-1,jj) ) * 0.5_wp230 zcmo(ji,jj,8) = zindb *( v_ice(ji,jj) * tmv(ji,jj) + v_ice(ji,jj-1) * tmv(ji,jj-1) ) * 0.5_wp230 zcmo(ji,jj,2) = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zinda 231 zcmo(ji,jj,3) = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zinda 232 zcmo(ji,jj,4) = diag_bot_gr(ji,jj) * 86400.0 ! Bottom thermodynamic ice production 233 zcmo(ji,jj,5) = diag_dyn_gr(ji,jj) * 86400.0 ! Dynamic ice production (rid/raft) 234 zcmo(ji,jj,22) = diag_lat_gr(ji,jj) * 86400.0 ! Lateral thermodynamic ice production 235 zcmo(ji,jj,23) = diag_sni_gr(ji,jj) * 86400.0 ! Snow ice production ice production 236 zcmo(ji,jj,24) = (tm_i(ji,jj) - rtt) * zinda 237 238 zcmo(ji,jj,6) = fbif(ji,jj)*at_i(ji,jj) 239 zcmo(ji,jj,7) = ( u_ice(ji,jj) * tmu(ji,jj) + u_ice(ji-1,jj) * tmu(ji-1,jj) ) * 0.5_wp 240 zcmo(ji,jj,8) = ( v_ice(ji,jj) * tmv(ji,jj) + v_ice(ji,jj-1) * tmv(ji,jj-1) ) * 0.5_wp 231 241 zcmo(ji,jj,9) = sst_m(ji,jj) 232 242 zcmo(ji,jj,10) = sss_m(ji,jj) … … 242 252 zcmo(ji,jj,19) = sprecip(ji,jj) 243 253 zcmo(ji,jj,20) = smt_i(ji,jj) 244 zcmo(ji,jj,21) = ot_i(ji,jj)245 254 zcmo(ji,jj,25) = et_i(ji,jj) 246 255 zcmo(ji,jj,26) = et_s(ji,jj) … … 249 258 250 259 zcmo(ji,jj,30) = bv_i(ji,jj) 251 zcmo(ji,jj,31) = hicol(ji,jj) 260 zcmo(ji,jj,31) = hicol(ji,jj) * zindb 252 261 zcmo(ji,jj,32) = strength(ji,jj) 253 262 zcmo(ji,jj,33) = SQRT( zcmo(ji,jj,7)*zcmo(ji,jj,7) + zcmo(ji,jj,8)*zcmo(ji,jj,8) ) 254 zcmo(ji,jj,34) = diag_sur_me(ji,jj) * 86400.0 * zinda! Surface melt255 zcmo(ji,jj,35) = diag_bot_me(ji,jj) * 86400.0 * zinda! Bottom melt263 zcmo(ji,jj,34) = diag_sur_me(ji,jj) * 86400.0 ! Surface melt 264 zcmo(ji,jj,35) = diag_bot_me(ji,jj) * 86400.0 ! Bottom melt 256 265 zcmo(ji,jj,36) = divu_i(ji,jj) 257 266 zcmo(ji,jj,37) = shear_i(ji,jj) 258 END DO 267 zcmo(ji,jj,38) = diag_res_pr(ji,jj) * 86400.0 ! Bottom melt 268 zcmo(ji,jj,39) = vt_i(ji,jj) ! ice volume 269 zcmo(ji,jj,40) = vt_s(ji,jj) ! snow volume 270 271 zcmo(ji,jj,41) = fsalt_rpo(ji,jj) 272 zcmo(ji,jj,42) = fsalt_res(ji,jj) 273 274 zcmo(ji,jj,43) = diag_trp_vi(ji,jj) * 86400.0 ! transport of ice volume 275 276 END DO 259 277 END DO 260 278 … … 285 303 ENDIF 286 304 305 CALL iom_put ('iceconc', zcmo(:,:,1) ) ! field1: ice concentration 306 CALL iom_put ('icethic_cea', zcmo(:,:,2) ) ! field2: ice thickness (i.e. icethi(:,:)) 307 CALL iom_put ('snowthic_cea', zcmo(:,:,3)) ! field3: snow thickness 308 CALL iom_put ('icebopr', zcmo(:,:,4) ) ! field4: daily bottom thermo ice production 309 CALL iom_put ('icedypr', zcmo(:,:,5) ) ! field5: daily dynamic ice production 310 CALL iom_put ('ioceflxb', zcmo(:,:,6) ) ! field6: Oceanic flux at the ice base 311 CALL iom_put ('uice_ipa', zcmo(:,:,7) ) ! field7: ice velocity u component 312 CALL iom_put ('vice_ipa', zcmo(:,:,8) ) ! field8: ice velocity v component 313 CALL iom_put ('isst', zcmo(:,:,9) ) ! field 9: sea surface temperature 314 CALL iom_put ('isss', zcmo(:,:,10) ) ! field 10: sea surface salinity 315 CALL iom_put ('qt_oc', zcmo(:,:,11) ) ! field 11: total flux at ocean surface 316 CALL iom_put ('qsr_oc', zcmo(:,:,12) ) ! field 12: solar flux at ocean surface 317 CALL iom_put ('qns_oc', zcmo(:,:,13) ) ! field 13: non-solar flux at ocean surface 318 !CALL iom_put ('hfbri', fhbri ) ! field 14: heat flux due to brine release 319 CALL iom_put ('qsr_io', zcmo(:,:,17) ) ! field 17: solar flux at ice/ocean surface 320 CALL iom_put ('qns_io', zcmo(:,:,18) ) ! field 18: non-solar flux at ice/ocean surface 321 !CALL iom_put ('snowpre', zcmo(:,:,19) * 86400) ! field 19 :snow precip 322 CALL iom_put ('micesalt', zcmo(:,:,20) ) ! field 20 :mean ice salinity 323 CALL iom_put ('miceage', zcmo(:,:,21) / 365) ! field 21: mean ice age 324 CALL iom_put ('icelapr',zcmo(:,:,22) ) ! field 22: daily lateral thermo ice prod. 325 CALL iom_put ('icesipr',zcmo(:,:,23) ) ! field 23: daily snowice ice prod. 326 CALL iom_put ('micet', zcmo(:,:,24) ) ! field 24: mean ice temperature 327 CALL iom_put ('icehc', zcmo(:,:,25) ) ! field 25: ice total heat content 328 CALL iom_put ('isnowhc', zcmo(:,:,26) ) ! field 26: snow total heat content 329 CALL iom_put ('icest', zcmo(:,:,27) ) ! field 27: ice surface temperature 330 CALL iom_put ('fsbri', zcmo(:,:,28) * 86400 ) ! field 28: brine salt flux 331 CALL iom_put ('fseqv', zcmo(:,:,29) * 86400 ) ! field 29: equivalent FW salt flux 332 CALL iom_put ('ibrinv', zcmo(:,:,30) *100 ) ! field 30: brine volume 333 CALL iom_put ('icecolf', zcmo(:,:,31) ) ! field 31: frazil ice collection thickness 334 CALL iom_put ('icestr', zcmo(:,:,32) * 0.001 ) ! field 32: ice strength 335 CALL iom_put ('icevel', zcmo(:,:,33) ) ! field 33: ice velocity 336 CALL iom_put ('isume', zcmo(:,:,34) ) ! field 34: surface melt 337 CALL iom_put ('ibome', zcmo(:,:,35) ) ! field 35: bottom melt 338 CALL iom_put ('idive', zcmo(:,:,36) * 1.0e8) ! field 36: divergence 339 CALL iom_put ('ishear', zcmo(:,:,37) * 1.0e8 ) ! field 37: shear 340 CALL iom_put ('icerepr', zcmo(:,:,38) ) ! field 38: daily prod./melting due to limupdate 341 CALL iom_put ('icevolu', zcmo(:,:,39) ) ! field 39: ice volume 342 CALL iom_put ('snowvol', zcmo(:,:,40) ) ! field 40: snow volume 343 CALL iom_put ('fsrpo', zcmo(:,:,41) * 86400 ) ! field 41: salt flux from ridging rafting 344 CALL iom_put ('fsres', zcmo(:,:,42) * 86400 ) ! field 42: salt flux from limupdate (resultant) 345 CALL iom_put ('icetrp', zcmo(:,:,43) ) ! field 43: ice volume transport 346 287 347 !----------------------------- 288 348 ! Thickness distribution file … … 302 362 DO jj = 1, jpj 303 363 DO ji = 1, jpi 304 zinda = MAX( zzero , SIGN( zone , a_i(ji,jj,jl) - 1.0e-6 ) )305 zoi(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , 1.0e-6 ) * zinda364 zinda = MAX( zzero , SIGN( zone , a_i(ji,jj,jl) - epsi06 ) ) 365 zoi(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi06 ) * zinda 306 366 END DO 307 367 END DO … … 314 374 DO jj = 1, jpj 315 375 DO ji = 1, jpi 316 zinda = MAX( zzero , SIGN( zone , a_i(ji,jj,jl) - 1.0e-6 ) )376 zinda = MAX( zzero , SIGN( zone , a_i(ji,jj,jl) - epsi06 ) ) 317 377 zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0* & 318 ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rtt ), - 1.0e-6 ) ) * &378 ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rtt ), - epsi06 ) ) * & 319 379 zinda / nlay_i 320 380 END DO … … 348 408 CALL wrk_dealloc( jpi, jpj, jpnoumax, zcmo, zcmoa ) 349 409 CALL wrk_dealloc( jpi, jpj, jpl, zmaskitd, zoi, zei ) 410 411 IF( nn_timing == 1 ) CALL timing_stop('limwri') 350 412 351 413 END SUBROUTINE lim_wri … … 381 443 field_25, field_26, field_27, field_28, field_29, field_30, & 382 444 field_31, field_32, field_33, field_34, field_35, field_36, & 383 field_37 445 field_37, field_38, field_39, field_40, field_41, field_42, field_43 384 446 385 447 TYPE(FIELD) , DIMENSION(jpnoumax) :: zfield … … 392 454 field_25, field_26, field_27, field_28, field_29, field_30, & 393 455 field_31, field_32, field_33, field_34, field_35, field_36, & 394 field_37, add_diag_swi456 field_37, field_38, field_39, field_40, field_41, field_42, field_43, add_diag_swi 395 457 !!------------------------------------------------------------------- 396 458 … … 435 497 zfield(36) = field_36 436 498 zfield(37) = field_37 499 zfield(38) = field_38 500 zfield(39) = field_39 501 zfield(40) = field_40 502 zfield(41) = field_41 503 zfield(42) = field_42 504 zfield(43) = field_43 437 505 438 506 DO nf = 1, noumef … … 460 528 ! 461 529 END SUBROUTINE lim_wri_init 530 531 SUBROUTINE lim_wri_state( kt, kid, kh_i ) 532 !!--------------------------------------------------------------------- 533 !! *** ROUTINE lim_wri_state *** 534 !! 535 !! ** Purpose : create a NetCDF file named cdfile_name which contains 536 !! the instantaneous ice state and forcing fields for ice model 537 !! Used to find errors in the initial state or save the last 538 !! ocean state in case of abnormal end of a simulation 539 !! 540 !! History : 541 !! 4.1 ! 2013-06 (C. Rousset) 542 !!---------------------------------------------------------------------- 543 INTEGER, INTENT( in ) :: kt ! ocean time-step index) 544 INTEGER, INTENT( in ) :: kid , kh_i 545 !!---------------------------------------------------------------------- 546 !CALL histvert( kid, "icethi", "L levels","m", jpl , hi_mean , nz ) 547 548 CALL histdef( kid, "iicethic", "Ice thickness" , "m" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 549 CALL histdef( kid, "iiceconc", "Ice concentration" , "%" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 550 CALL histdef( kid, "iicetemp", "Ice temperature" , "C" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 551 CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 552 CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 553 CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 554 CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 555 CALL histdef( kid, "iicesflx", "Solar flux over ocean" , "w/m2" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 556 CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 557 CALL histdef( kid, "isnowpre", "Snow precipitation" , "kg/m2/s", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 558 CALL histdef( kid, "iicesali", "Ice salinity" , "PSU" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 559 CALL histdef( kid, "iicevolu", "Ice volume" , "m" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 560 CALL histdef( kid, "iicedive", "Ice divergence" , "10-8s-1", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 561 562 !CALL histdef( kid, "iice_itd", "Ice concentration by cat", "%" , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 563 !CALL histdef( kid, "iice_hid", "Ice thickness by cat" , "m" , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 564 !CALL histdef( kid, "iice_hsd", "Snow thickness by cat" , "m" , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 565 !CALL histdef( kid, "iice_std", "Ice salinity by cat" , "PSU" , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 566 567 CALL histend( kid, snc4set ) ! end of the file definition 568 569 CALL histwrite( kid, "iicethic", kt, icethi , jpi*jpj, (/1/) ) 570 CALL histwrite( kid, "iiceconc", kt, at_i , jpi*jpj, (/1/) ) 571 CALL histwrite( kid, "iicetemp", kt, tm_i - rtt , jpi*jpj, (/1/) ) 572 CALL histwrite( kid, "iicevelu", kt, u_ice , jpi*jpj, (/1/) ) 573 CALL histwrite( kid, "iicevelv", kt, v_ice , jpi*jpj, (/1/) ) 574 CALL histwrite( kid, "iicestru", kt, utau_ice , jpi*jpj, (/1/) ) 575 CALL histwrite( kid, "iicestrv", kt, vtau_ice , jpi*jpj, (/1/) ) 576 CALL histwrite( kid, "iicesflx", kt, qsr , jpi*jpj, (/1/) ) 577 CALL histwrite( kid, "iicenflx", kt, qns , jpi*jpj, (/1/) ) 578 CALL histwrite( kid, "isnowpre", kt, sprecip , jpi*jpj, (/1/) ) 579 CALL histwrite( kid, "iicesali", kt, smt_i , jpi*jpj, (/1/) ) 580 CALL histwrite( kid, "iicevolu", kt, vt_i , jpi*jpj, (/1/) ) 581 CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8 , jpi*jpj, (/1/) ) 582 583 !CALL histwrite( kid, "iice_itd", kt, a_i , jpi*jpj*jpl, (/1/) ) ! area 584 !CALL histwrite( kid, "iice_hid", kt, ht_i , jpi*jpj*jpl, (/1/) ) ! thickness 585 !CALL histwrite( kid, "iice_hsd", kt, ht_s , jpi*jpj*jpl, (/1/) ) ! snow depth 586 !CALL histwrite( kid, "iice_std", kt, sm_i , jpi*jpj*jpl, (/1/) ) ! salinity 587 588 END SUBROUTINE lim_wri_state 462 589 463 590 #else -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90
r2715 r3938 111 111 zcmo(ji,jj,13) = qns(ji,jj) 112 112 ! See thersf for the coefficient 113 zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce113 zcmo(ji,jj,14) = - emps(ji,jj) * rday ! converted in Kg/m2/day = mm/day 114 114 zcmo(ji,jj,15) = utau_ice(ji,jj) 115 115 zcmo(ji,jj,16) = vtau_ice(ji,jj) … … 154 154 rcmoy(ji,jj,13) = qns(ji,jj) 155 155 ! See thersf for the coefficient 156 rcmoy(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce156 rcmoy(ji,jj,14) = - emps(ji,jj) * rday ! converted in Kg/m2/day = mm/day 157 157 rcmoy(ji,jj,15) = utau_ice(ji,jj) 158 158 rcmoy(ji,jj,16) = vtau_ice(ji,jj) -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r2715 r3938 22 22 REAL(wp), PUBLIC :: hicmin = 0.2 !: (REMOVE) 23 23 REAL(wp), PUBLIC :: hiclim = 0.05 !: minimum ice thickness 24 REAL(wp), PUBLIC :: amax = 0.999 !: maximum lead fraction25 24 REAL(wp), PUBLIC :: sbeta = 1.0 !: numerical scheme for diffusion in ice (REMOVE) 26 25 REAL(wp), PUBLIC :: parlat = 0.0 !: (REMOVE) … … 83 82 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: i0 !: fraction of radiation transmitted to the ice 84 83 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: old_ht_i_b !: Ice thickness at the beginnning of the time step [m] 85 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::old_ht_s_b !: Snow thickness at the beginning of the time step [m]84 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: old_ht_s_b !: Snow thickness at the beginning of the time step [m] 86 85 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fsbri_1d !: Salt flux due to brine drainage 87 86 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhbri_1d !: Heat flux due to brine drainage … … 108 107 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: o_i_b !: Ice age [days] 109 108 109 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: iatte_1d !: clem attenuation coef of the input solar flux (unitless) 110 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: oatte_1d !: clem attenuation coef of the input solar flux (unitless) 111 110 112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_b !: corresponding to the 2D var t_s 111 113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_i_b !: corresponding to the 2D var t_i … … 157 159 & fltbif_1d(jpij) , fscbq_1d (jpij) , qsr_ice_1d (jpij) , & 158 160 & fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qnsr_ice_1d(jpij) , & 159 & qfvbq_1d (jpij) , t_bo_b (jpij) , STAT=ierr(1) ) 161 & qfvbq_1d (jpij) , t_bo_b (jpij) , iatte_1d (jpij) , & 162 & oatte_1d (jpij) , STAT=ierr(1) ) 160 163 ! 161 164 ALLOCATE( sprecip_1d (jpij) , frld_1d (jpij) , at_i_b (jpij) , & -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
r3294 r3938 8 8 !! 3.3 ! 2010-09 (D. Storkey) add ice boundary conditions 9 9 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge 10 !! - ! 2012-01 (C. Rousset) add ice boundary conditions for lim3 10 11 !!---------------------------------------------------------------------- 11 12 #if defined key_bdy … … 44 45 REAL, POINTER, DIMENSION(:) :: hicif 45 46 REAL, POINTER, DIMENSION(:) :: hsnif 47 #elif defined key_lim3 48 REAL, POINTER, DIMENSION(:,:) :: a_i !: now ice leads fraction climatology 49 REAL, POINTER, DIMENSION(:,:) :: ht_i !: Now ice thickness climatology 50 REAL, POINTER, DIMENSION(:,:) :: ht_s !: now snow thickness 46 51 #endif 47 52 END TYPE OBC_DATA … … 73 78 INTEGER, DIMENSION(jp_bdy) :: nn_tra_dta !: = 0 use the initial state as bdy dta ; 74 79 !: = 1 read it in a NetCDF file 75 #if defined key_lim276 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim 2! Choice of boundary condition for sea ice variables77 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim 2_dta!: = 0 use the initial state as bdy dta ;80 #if ( defined key_lim2 || defined key_lim3 ) 81 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim ! Choice of boundary condition for sea ice variables 82 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim_dta !: = 0 use the initial state as bdy dta ; 78 83 !: = 1 read it in a NetCDF file 79 84 #endif -
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r3294 r3938 11 11 !! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions 12 12 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge 13 !! - ! 2012-01 (C. Rousset) add ice boundary conditions for lim3 13 14 !!---------------------------------------------------------------------- 14 15 #if defined key_bdy … … 31 32 #if defined key_lim2 32 33 USE ice_2 34 #elif defined key_lim3 35 USE par_ice 36 USE ice 37 USE limcat_1D ! redistribute ice input into categories 33 38 #endif 34 39 … … 48 53 49 54 TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr ! array of pointers to nbmap 55 56 #if defined key_lim3 57 LOGICAL :: ll_bdylim3 ! determine whether ice input is lim2 (F) or lim3 (T) type