Changeset 4045
- Timestamp:
- 2013-09-25T16:38:24+02:00 (11 years ago)
- Location:
- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO
- Files:
-
- 5 added
- 29 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90
r3625 r4045 31 31 32 32 !!---------------------------------------------------------------------- 33 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)33 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 34 34 !! $Id$ 35 35 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r3625 r4045 188 188 REAL(wp), PUBLIC :: alphaevp = 1._wp !: coeficient of the internal stresses !SB 189 189 REAL(wp), PUBLIC :: unit_fac = 1.e+09_wp !: conversion factor for ice / snow enthalpy 190 REAL(wp), PUBLIC :: hminrhg = 0.05_wp !: clem : ice thickness (in m) below which ice velocity is set to ocean velocity 190 191 191 192 ! !!** ice-salinity namelist (namicesal) ** … … 395 396 LOGICAL , PUBLIC :: ln_limdyn = .TRUE. !: flag for ice dynamics (T) or not (F) 396 397 LOGICAL , PUBLIC :: ln_nicep = .TRUE. !: flag for sea-ice points output (T) or not (F) 397 REAL(wp) , PUBLIC :: hsndif = 0.e0 !: computation of temp. in snow (0) or not (9999)398 REAL(wp) , PUBLIC :: hicdif = 0.e0 !: computation of temp. in ice (0) or not (9999)399 398 REAL(wp) , PUBLIC :: cai = 1.40e-3 !: atmospheric drag over sea ice 400 399 REAL(wp) , PUBLIC :: cao = 1.00e-3 !: atmospheric drag over ocean 401 REAL(wp) , DIMENSION(2), PUBLIC :: acrit = (/ 1.e-06 , 1.e-06 /) !: minimum fraction for leads in402 ! ! : north and south hemisphere400 REAL(wp) , PUBLIC :: amax = 0.99 !: maximum ice concentration 401 ! ! 403 402 !!-------------------------------------------------------------------------- 404 403 !! * Ice diagnostics 405 404 !!-------------------------------------------------------------------------- 406 405 !! Check if everything down here is necessary 406 LOGICAL , PUBLIC :: ln_limdiahsb = .TRUE. !: flag for ice diag (T) or not (F) 407 407 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_newice !: volume of ice formed in the leads 408 408 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dv_dt_thd !: thermodynamic growth rates … … 414 414 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_bot_me ! vertical bottom melt 415 415 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_sur_me ! vertical surface melt 416 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_res_pr ! production (growth+melt) due to limupdate 417 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vi ! transport of ice volume 416 418 INTEGER , PUBLIC :: jiindx, jjindx !: indexes of the debugging point 417 419 418 420 !!---------------------------------------------------------------------- 419 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2010)421 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 420 422 !! $Id$ 421 423 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 527 529 & izero (jpi,jpj,jpl) , diag_bot_gr(jpi,jpj) , diag_dyn_gr(jpi,jpj) , & 528 530 & fstroc (jpi,jpj,jpl) , diag_bot_me(jpi,jpj) , diag_sur_me(jpi,jpj) , & 529 & fhbricat (jpi,jpj,jpl) , v_newice (jpi,jpj), STAT=ierr(ii) )531 & fhbricat (jpi,jpj,jpl) , diag_res_pr(jpi,jpj) , diag_trp_vi(jpi,jpj) , v_newice(jpi,jpj) , STAT=ierr(ii) ) 530 532 531 533 ice_alloc = MAXVAL( ierr(:) ) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
r3625 r4045 40 40 41 41 !!---------------------------------------------------------------------- 42 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)42 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 43 43 !! $Id$ 44 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 126 126 !! ** input : Namelist namicerun 127 127 !!------------------------------------------------------------------- 128 NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, a crit, hsndif, hicdif, cai, cao, ln_nicep128 NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb 129 129 !!------------------------------------------------------------------- 130 130 ! … … 142 142 WRITE(numout,*) ' ~~~~~~' 143 143 WRITE(numout,*) ' switch for ice dynamics (1) or not (0) ln_limdyn = ', ln_limdyn 144 WRITE(numout,*) ' minimum fraction for leads in the NH (SH) acrit(1/2) = ', acrit(:) 145 WRITE(numout,*) ' computation of temp. in snow (=0) or not (=9999) hsndif = ', hsndif 146 WRITE(numout,*) ' computation of temp. in ice (=0) or not (=9999) hicdif = ', hicdif 144 WRITE(numout,*) ' maximum ice concentration = ', amax 147 145 WRITE(numout,*) ' atmospheric drag over sea ice = ', cai 148 146 WRITE(numout,*) ' atmospheric drag over ocean = ', cao 149 147 WRITE(numout,*) ' Several ice points in the ice or not in ocean.output = ', ln_nicep 148 WRITE(numout,*) ' Diagnose heat/salt budget or not ln_limdiahsb = ', ln_limdiahsb 150 149 ENDIF 151 150 ! -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
r3764 r4045 15 15 !! lim_adv_y : advection of sea ice on y axis 16 16 !!---------------------------------------------------------------------- 17 USE dom_oce ! ocean domain 18 USE ice ! LIM-3 variables 19 USE dom_ice ! LIM-3 domain 20 USE lbclnk ! lateral boundary condition - MPP exchanges 21 USE in_out_manager ! I/O manager 22 USE prtctl ! Print control 23 USE lib_mpp ! MPP library 24 USE wrk_nemo ! work arrays 25 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 26 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 17 USE dom_oce ! ocean domain 18 USE dom_ice ! LIM-3 domain 19 USE ice ! LIM-3 variables 20 USE lbclnk ! lateral boundary condition - MPP exchanges 21 USE in_out_manager ! I/O manager 22 USE prtctl ! Print control 23 USE lib_mpp ! MPP library 24 USE wrk_nemo ! work arrays 25 USE lib_fortran ! to use key_nosignedzero 27 26 28 27 IMPLICIT NONE … … 39 38 # include "vectopt_loop_substitute.h90" 40 39 !!---------------------------------------------------------------------- 41 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)40 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 42 41 !! $Id$ 43 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r3625 r4045 30 30 31 31 !!---------------------------------------------------------------------- 32 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)32 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 33 33 !! $Id$ 34 34 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r3625 r4045 15 15 !! lim_dyn_init : initialization and namelist read 16 16 !!---------------------------------------------------------------------- 17 USE phycst ! physical constants18 USE dom_oce ! ocean space and time domain19 USE sbc_oce ! Surface boundary condition: ocean fields20 USE sbc_ice ! Surface boundary condition: ice fields21 USE ice ! LIM-3 variables22 USE par_ice ! LIM-3 parameters23 USE dom_ice ! LIM-3 domain24 USE limrhg ! LIM-3 rheology25 USE lbclnk ! lateral boundary conditions - MPP exchanges26 USE lib_mpp ! MPP library27 USE wrk_nemo ! work arrays28 USE in_out_manager ! I/O manager29 USE prtctl ! Print control30 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)17 USE phycst ! physical constants 18 USE dom_oce ! ocean space and time domain 19 USE sbc_oce ! Surface boundary condition: ocean fields 20 USE sbc_ice ! Surface boundary condition: ice fields 21 USE ice ! LIM-3 variables 22 USE par_ice ! LIM-3 parameters 23 USE dom_ice ! LIM-3 domain 24 USE limrhg ! LIM-3 rheology 25 USE lbclnk ! lateral boundary conditions - MPP exchanges 26 USE lib_mpp ! MPP library 27 USE wrk_nemo ! work arrays 28 USE in_out_manager ! I/O manager 29 USE prtctl ! Print control 30 USE lib_fortran ! glob_sum 31 31 32 32 IMPLICIT NONE … … 38 38 # include "vectopt_loop_substitute.h90" 39 39 !!---------------------------------------------------------------------- 40 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)40 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 41 41 !! $Id$ 42 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 65 65 REAL(wp), POINTER, DIMENSION(:) :: zmsk ! i-averaged of tmask 66 66 REAL(wp), POINTER, DIMENSION(:,:) :: zu_io, zv_io ! ice-ocean velocity 67 !!--------------------------------------------------------------------- 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 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 69 !!--------------------------------------------------------------------- 68 70 69 71 CALL wrk_alloc( jpi, jpj, zu_io, zv_io ) 70 72 CALL wrk_alloc( jpj, zind, zmsk ) 73 74 ! ------------------------------- 75 !- check conservation (C Rousset) 76 IF (ln_limdiahsb) THEN 77 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 78 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 79 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 80 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 81 ENDIF 82 !- check conservation (C Rousset) 83 ! ------------------------------- 71 84 72 85 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) … … 208 221 ENDIF 209 222 ! 223 224 ! ------------------------------- 225 !- check conservation (C Rousset) 226 IF (ln_limdiahsb) THEN 227 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 228 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 229 230 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 231 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 232 233 zchk_vmin = glob_min(v_i) 234 zchk_amax = glob_max(SUM(a_i,dim=3)) 235 zchk_amin = glob_min(a_i) 236 237 IF(lwp) THEN 238 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limdyn) = ',(zchk_v_i * rday) 239 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * rday) 240 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limdyn) = ',(zchk_vmin * 1.e-3) 241 !IF ( zchk_amax > amax+1.e-10 ) WRITE(numout,*) 'violation a_i>amax (limdyn) = ',zchk_amax 242 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limdyn) = ',zchk_amin 243 ENDIF 244 ENDIF 245 !- check conservation (C Rousset) 246 ! ------------------------------- 247 210 248 CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 211 249 CALL wrk_dealloc( jpj, zind, zmsk ) … … 229 267 & dm, nbiter, nbitdr, om, resl, cw, angvg, pstar, & 230 268 & c_rhg, etamn, creepl, ecc, ahi0, & 231 & nevp, telast, alphaevp 269 & nevp, telast, alphaevp, hminrhg 232 270 !!------------------------------------------------------------------- 233 271 … … 257 295 WRITE(numout,*) ' timescale for elastic waves telast = ', telast 258 296 WRITE(numout,*) ' coefficient for the solution of int. stresses alphaevp = ', alphaevp 297 WRITE(numout,*) ' min ice thickness for rheology calculations hminrhg = ', hminrhg 259 298 ENDIF 260 299 ! -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r3625 r4045 35 35 # include "vectopt_loop_substitute.h90" 36 36 !!---------------------------------------------------------------------- 37 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2010)37 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 38 38 !! $Id$ 39 39 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 137 137 END DO ! end of sub-time step loop 138 138 139 ! ----------------------- 140 !!! final step (clem) !!! 141 DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction 142 DO ji = 1 , fs_jpim1 ! vector opt. 143 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 144 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 145 END DO 146 END DO 147 ! 148 DO jj= 2, jpjm1 ! diffusive trend : divergence of the fluxes 149 DO ji = fs_2 , fs_jpim1 ! vector opt. 150 zdiv (ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj ) & 151 & + zflv(ji,jj) - zflv(ji ,jj-1) ) / ( e1t (ji,jj) * e2t (ji,jj) ) 152 ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) ) 153 END DO 154 END DO 155 CALL lbc_lnk( ptab, 'T', 1. ) ! lateral boundary condition 156 !!! final step (clem) !!! 157 ! ----------------------- 158 139 159 IF(ln_ctl) THEN 140 160 zrlx(:,:) = ptab(:,:) - ztab0(:,:) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r3625 r4045 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) 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 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r3764 r4045 5 5 !!====================================================================== 6 6 !! History : LIM ! 2006-02 (M. Vancoppenolle) Original code 7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & fsalt_rpo7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_mec 8 8 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 9 9 !!---------------------------------------------------------------------- … … 12 12 !! 'key_lim3' LIM-3 sea-ice model 13 13 !!---------------------------------------------------------------------- 14 USE par_oce ! ocean parameters 15 USE dom_oce ! ocean domain 16 USE phycst ! physical constants (ocean directory) 17 USE sbc_oce ! surface boundary condition: ocean fields 18 USE thd_ice ! LIM thermodynamics 19 USE ice ! LIM variables 20 USE par_ice ! LIM parameters 21 USE dom_ice ! LIM domain 22 USE limthd_lac ! LIM 23 USE limvar ! LIM 24 USE limcons ! LIM 25 USE in_out_manager ! I/O manager 26 USE lbclnk ! lateral boundary condition - MPP exchanges 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! work arrays 29 USE prtctl ! Print control 30 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 14 USE par_oce ! ocean parameters 15 USE dom_oce ! ocean domain 16 USE phycst ! physical constants (ocean directory) 17 USE sbc_oce ! surface boundary condition: ocean fields 18 USE thd_ice ! LIM thermodynamics 19 USE ice ! LIM variables 20 USE par_ice ! LIM parameters 21 USE dom_ice ! LIM domain 22 USE limthd_lac ! LIM 23 USE limvar ! LIM 24 USE limcons ! LIM 25 USE in_out_manager ! I/O manager 26 USE lbclnk ! lateral boundary condition - MPP exchanges 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! work arrays 29 USE prtctl ! Print control 30 ! Check budget (Rousset) 31 USE iom ! I/O manager 32 USE lib_fortran ! glob_sum 33 USE limdiahsb 32 34 33 35 IMPLICIT NONE … … 62 64 REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier 63 65 REAL(wp), PARAMETER :: kraft = 2.0_wp ! rafting multipliyer 66 REAL(wp), PARAMETER :: kamax = 1.0 64 67 65 68 REAL(wp) :: Cp ! … … 141 144 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2) 142 145 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 146 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) 147 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 148 ! mass and salt flux (clem) 149 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... 143 150 !!----------------------------------------------------------------------------- 144 151 145 152 CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 153 154 CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem 146 155 147 156 IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only) … … 151 160 CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ') 152 161 ENDIF 162 163 IF( ln_limdyn ) THEN ! Start ridging and rafting ! 164 ! ------------------------------- 165 !- check conservation (C Rousset) 166 IF (ln_limdiahsb) THEN 167 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 168 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 169 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 170 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 171 ENDIF 172 !- check conservation (C Rousset) 173 ! ------------------------------- 174 175 ! mass and salt flux init (clem) 176 zviold(:,:,:) = v_i(:,:,:) 177 zvsold(:,:,:) = v_s(:,:,:) 178 zsmvold(:,:,:) = smv_i(:,:,:) 153 179 154 180 !-----------------------------------------------------------------------------! … … 204 230 ! to give asum = 1.0 after ridging. 205 231 206 divu_adv(ji,jj) = ( 1._wp- asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep232 divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep 207 233 208 234 IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) … … 288 314 DO jj = 1, jpj 289 315 DO ji = 1, jpi 290 IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN316 IF (ABS(asum(ji,jj) - kamax ) .LT. epsi11) THEN 291 317 closing_net(ji,jj) = 0._wp 292 318 opning (ji,jj) = 0._wp 293 319 ELSE 294 320 iterate_ridging = 1 295 divu_adv (ji,jj) = ( 1._wp - asum(ji,jj)) * r1_rdtice321 divu_adv (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice 296 322 closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 297 323 opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) ) … … 330 356 DO ji = 1, jpi 331 357 332 IF( ABS( asum(ji,jj) - 1.0 ) > epsi11 )asum_error = .true.358 IF(ABS(asum(ji,jj) - kamax) > epsi11 ) asum_error = .true. 333 359 334 360 dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice … … 349 375 DO jj = 1, jpj 350 376 DO ji = 1, jpi 351 IF( ABS( asum(ji,jj) - 1._wp) > epsi11 ) THEN ! there is a bug377 IF( ABS( asum(ji,jj) - kamax) > epsi11 ) THEN ! there is a bug 352 378 WRITE(numout,*) ' ' 353 379 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) … … 377 403 CALL lim_var_glo2eqv 378 404 CALL lim_itd_me_zapsmall 405 406 !-------------------------------- 407 ! Update mass/salt fluxes (clem) 408 !-------------------------------- 409 DO jl = 1, jpl 410 DO jj = 1, jpj 411 DO ji = 1, jpi 412 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 413 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 414 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 415 sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice 416 END DO 417 END DO 418 END DO 379 419 380 420 !----------------- … … 425 465 ENDIF 426 466 467 ! ------------------------------- 468 !- check conservation (C Rousset) 469 IF (ln_limdiahsb) THEN 470 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 471 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 472 473 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 474 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 475 476 zchk_vmin = glob_min(v_i) 477 zchk_amax = glob_max(SUM(a_i,dim=3)) 478 zchk_amin = glob_min(a_i) 479 480 IF(lwp) THEN 481 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_me) = ',(zchk_v_i * rday) 482 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday) 483 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_me) = ',(zchk_vmin * 1.e-3) 484 IF ( zchk_amax > kamax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_me) = ',zchk_amax 485 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_me) = ',zchk_amin 486 ENDIF 487 ENDIF 488 !- check conservation (C Rousset) 489 ! ------------------------------- 490 427 491 !-------------------------! 428 492 ! Back to initial values … … 448 512 449 513 ! 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 ) 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 514 !clem@order 515 ! DO jl = 1, jpl 516 ! DO jk = 1, nlay_i 517 ! DO jj = 1, jpj 518 ! DO ji = 1, jpi 519 ! IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 520 ! ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 521 ! old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl) 522 ! d_e_i_trp(ji,jj,jk,jl) = 0._wp 523 ! ENDIF 524 ! END DO 525 ! END DO 526 ! END DO 527 ! END DO 528 ! 529 ! DO jl = 1, jpl 530 ! DO jj = 1, jpj 531 ! DO ji = 1, jpi 532 ! IF( old_v_i (ji,jj,jl) < epsi06 .AND. & 533 ! d_v_i_trp(ji,jj,jl) > epsi06 ) THEN 534 ! old_v_i (ji,jj,jl) = d_v_i_trp(ji,jj,jl) 535 ! d_v_i_trp (ji,jj,jl) = 0._wp 536 ! old_a_i (ji,jj,jl) = d_a_i_trp(ji,jj,jl) 537 ! d_a_i_trp (ji,jj,jl) = 0._wp 538 ! old_v_s (ji,jj,jl) = d_v_s_trp(ji,jj,jl) 539 ! d_v_s_trp (ji,jj,jl) = 0._wp 540 ! old_e_s (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 541 ! d_e_s_trp (ji,jj,1,jl) = 0._wp 542 ! old_oa_i (ji,jj,jl) = d_oa_i_trp(ji,jj,jl) 543 ! d_oa_i_trp(ji,jj,jl) = 0._wp 544 ! IF( num_sal == 2 ) old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 545 ! d_smv_i_trp(ji,jj,jl) = 0._wp 546 ! ENDIF 547 ! END DO 548 ! END DO 549 ! END DO 550 !clem@order 551 ENDIF ! ln_limdyn=.true. 485 552 ! 486 553 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 554 ! 555 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem 487 556 ! 488 557 END SUBROUTINE lim_itd_me … … 1086 1155 afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 1087 1156 1088 IF (afrac(ji,jj) > 1.0+ epsi11) THEN !riging1157 IF (afrac(ji,jj) > kamax + epsi11) THEN !riging 1089 1158 large_afrac = .true. 1090 ELSEIF (afrac(ji,jj) > 1.0) THEN ! roundoff error1091 afrac(ji,jj) = 1.01159 ELSEIF (afrac(ji,jj) > kamax) THEN ! roundoff error 1160 afrac(ji,jj) = kamax 1092 1161 ENDIF 1093 IF (afrft(ji,jj) > 1.0+ epsi11) THEN !rafting1162 IF (afrft(ji,jj) > kamax + epsi11) THEN !rafting 1094 1163 large_afrft = .true. 1095 ELSEIF (afrft(ji,jj) > 1.0) THEN ! roundoff error1096 afrft(ji,jj) = 1.01164 ELSEIF (afrft(ji,jj) > kamax) THEN ! roundoff error 1165 afrft(ji,jj) = kamax 1097 1166 ENDIF 1098 1167 … … 1137 1206 1138 1207 ! ! excess of salt is flushed into the ocean 1139 sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice1140 1141 rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic / rau0 ! increase in ice volume du to seawater frozen in voids1142 1208 !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 1209 1210 !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic ! gurvan: increase in ice volume du to seawater frozen in voids 1211 1143 1212 !------------------------------------ 1144 1213 ! 3.6 Increment ridging diagnostics … … 1150 1219 dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 1151 1220 dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 1152 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice1221 !clem diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 1153 1222 opening (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 1154 1223 … … 1217 1286 1218 1287 ! 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_i1288 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 1220 1289 1221 1290 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) … … 1240 1309 ji = indxi(ij) 1241 1310 jj = indxj(ij) 1242 IF( afrac(ji,jj) > 1.0+ epsi11 ) THEN1311 IF( afrac(ji,jj) > kamax + epsi11 ) THEN 1243 1312 WRITE(numout,*) '' 1244 1313 WRITE(numout,*) ' ardg > a_i' … … 1252 1321 ji = indxi(ij) 1253 1322 jj = indxj(ij) 1254 IF( afrft(ji,jj) > 1.0+ epsi11 ) THEN1323 IF( afrft(ji,jj) > kamax + epsi11 ) THEN 1255 1324 WRITE(numout,*) '' 1256 1325 WRITE(numout,*) ' arft > a_i' -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r3764 r4045 19 19 !! lim_itd_shiftice : 20 20 !!---------------------------------------------------------------------- 21 USE par_oce ! ocean parameters 22 USE dom_oce ! ocean domain 23 USE phycst ! physical constants (ocean directory) 24 USE ice ! LIM-3 variables 25 USE par_ice ! LIM-3 parameters 26 USE dom_ice ! LIM-3 domain 27 USE thd_ice ! LIM-3 thermodynamic variables 28 USE limthd_lac ! LIM-3 lateral accretion 29 USE limvar ! LIM-3 variables 30 USE limcons ! LIM-3 conservation 31 USE prtctl ! Print control 32 USE in_out_manager ! I/O manager 33 USE lib_mpp ! MPP library 34 USE wrk_nemo ! work arrays 35 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 36 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 21 USE dom_ice ! LIM-3 domain 22 USE par_oce ! ocean parameters 23 USE dom_oce ! ocean domain 24 USE phycst ! physical constants (ocean directory) 25 USE thd_ice ! LIM-3 thermodynamic variables 26 USE ice ! LIM-3 variables 27 USE par_ice ! LIM-3 parameters 28 USE limthd_lac ! LIM-3 lateral accretion 29 USE limvar ! LIM-3 variables 30 USE limcons ! LIM-3 conservation 31 USE prtctl ! Print control 32 USE in_out_manager ! I/O manager 33 USE lib_mpp ! MPP library 34 USE wrk_nemo ! work arrays 35 USE lib_fortran ! to use key_nosignedzero 37 36 38 37 IMPLICIT NONE … … 45 44 PUBLIC lim_itd_shiftice 46 45 47 REAL(wp) :: epsi20 = 1 e-20_wp ! constant values48 REAL(wp) :: epsi13 = 1 e-13_wp !49 REAL(wp) :: epsi10 = 1 e-10_wp !46 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values 47 REAL(wp) :: epsi13 = 1.e-13_wp ! 48 REAL(wp) :: epsi10 = 1.e-10_wp ! 50 49 51 50 !!---------------------------------------------------------------------- 52 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2010)51 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 53 52 !! $Id$ 54 53 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 67 66 ! 68 67 INTEGER :: jl, ja, jm, jbnd1, jbnd2 ! ice types dummy loop index 69 70 !!------------------------------------------------------------------ 68 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) 69 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (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( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 77 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 78 ENDIF 79 !- check conservation (C Rousset) 80 ! ------------------------------- 71 81 72 82 IF( kt == nit000 .AND. lwp ) THEN … … 107 117 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 108 118 d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 109 119 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 110 120 d_smv_i_thd(:,:,:) = 0._wp 111 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 121 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 122 123 ! diag only (clem) 124 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 112 125 113 126 IF(ln_ctl) THEN ! Control print … … 142 155 END DO 143 156 ENDIF 144 157 ! 158 ! ------------------------------- 159 !- check conservation (C Rousset) 160 IF (ln_limdiahsb) THEN 161 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 162 zchk_fw = glob_sum( rdm_ice(:,:) * 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 ) ) * r1_rdtice 165 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 166 167 zchk_vmin = glob_min(v_i) 168 zchk_amax = glob_max(SUM(a_i,dim=3)) 169 zchk_amin = glob_min(a_i) 170 171 IF(lwp) THEN 172 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_th) = ',(zchk_v_i * rday) 173 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * rday) 174 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_th) = ',(zchk_vmin * 1.e-3) 175 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_th) = ',zchk_amax 176 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_th) = ',zchk_amin 177 ENDIF 178 ENDIF 179 !- check conservation (C Rousset) 180 ! ------------------------------- 181 ! 145 182 !- Recover Old values 146 a_i(:,:,:) = old_a_i (:,:,:)147 v_s(:,:,:) = old_v_s (:,:,:)148 v_i(:,:,:) = old_v_i (:,:,:)149 e_s(:,:,:,:) = old_e_s (:,:,:,:)150 e_i(:,:,:,:) = old_e_i (:,:,:,:)151 ! 183 a_i(:,:,:) = old_a_i (:,:,:) 184 v_s(:,:,:) = old_v_s (:,:,:) 185 v_i(:,:,:) = old_v_i (:,:,:) 186 e_s(:,:,:,:) = old_e_s (:,:,:,:) 187 e_i(:,:,:,:) = old_e_i (:,:,:,:) 188 !?? oa_i(:,:,:) = old_oa_i(:,:,:) 152 189 IF( num_sal == 2 ) smv_i(:,:,:) = old_smv_i(:,:,:) 153 190 ! … … 172 209 ! 173 210 INTEGER :: ji, jj, jl ! dummy loop index 174 INTEGER :: zji, zjj, nd ! local integer 211 INTEGER :: ii, ij ! 2D corresponding indices to ji 212 INTEGER :: nd ! local integer 175 213 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 176 REAL(wp) :: zx2, zwk2, zda0, zetamax , zhimin! - -214 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 177 215 REAL(wp) :: zx3, zareamin, zindb ! - - 178 216 CHARACTER (len = 15) :: fieldid … … 210 248 CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 211 249 212 zhimin = 0.1 !minimum ice thickness tolerated by the model213 250 zareamin = epsi10 !minimum area in thickness categories tolerated by the conceptors of the model 214 251 … … 240 277 DO jj = 1, jpj 241 278 DO ji = 1, jpi 242 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl) )) !0 if no ice and 1 if yes279 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 243 280 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi10) * zindb 244 zindb = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl) )) !0 if no ice and 1 if yes281 zindb = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 245 282 zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),epsi10) * zindb 246 283 IF( a_i(ji,jj,jl) > 1e-6 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) … … 285 322 DO jl = klbnd, kubnd - 1 286 323 DO ji = 1, nbrem 287 zji = nind_i(ji)288 zjj = nind_j(ji)324 ii = nind_i(ji) 325 ij = nind_j(ji) 289 326 ! 290 IF ( ( zht_i_o( zji,zjj,jl) .GT.epsi10 ) .AND. &291 ( zht_i_o( zji,zjj,jl+1).GT.epsi10 ) ) THEN327 IF ( ( zht_i_o(ii,ij,jl) .GT.epsi10 ) .AND. & 328 ( zht_i_o(ii,ij,jl+1).GT.epsi10 ) ) THEN 292 329 !interpolate between adjacent category growth rates 293 zslope = ( zdhice( zji,zjj,jl+1) - zdhice(zji,zjj,jl) ) / &294 ( zht_i_o ( zji,zjj,jl+1) - zht_i_o (zji,zjj,jl) )295 zhbnew( zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) + &296 zslope * ( hi_max(jl) - zht_i_o( zji,zjj,jl) )297 ELSEIF (zht_i_o( zji,zjj,jl).gt.epsi10) THEN298 zhbnew( zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl)299 ELSEIF (zht_i_o( zji,zjj,jl+1).gt.epsi10) THEN300 zhbnew( zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl+1)330 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / & 331 ( zht_i_o (ii,ij,jl+1) - zht_i_o (ii,ij,jl) ) 332 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + & 333 zslope * ( hi_max(jl) - zht_i_o(ii,ij,jl) ) 334 ELSEIF (zht_i_o(ii,ij,jl).gt.epsi10) THEN 335 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 336 ELSEIF (zht_i_o(ii,ij,jl+1).gt.epsi10) THEN 337 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 301 338 ELSE 302 zhbnew( zji,zjj,jl) = hi_max(jl)339 zhbnew(ii,ij,jl) = hi_max(jl) 303 340 ENDIF 304 341 END DO … … 307 344 DO ji = 1, nbrem 308 345 ! jl, ji 309 zji = nind_i(ji)310 zjj = nind_j(ji)346 ii = nind_i(ji) 347 ij = nind_j(ji) 311 348 ! jl, ji 312 IF ( ( a_i( zji,zjj,jl) .GT.epsi10) .AND. &313 ( ht_i( zji,zjj,jl).GE. zhbnew(zji,zjj,jl) ) &349 IF ( ( a_i(ii,ij,jl) .GT.epsi10) .AND. & 350 ( ht_i(ii,ij,jl).GE. zhbnew(ii,ij,jl) ) & 314 351 ) THEN 315 zremap_flag( zji,zjj) = 0316 ELSEIF ( ( a_i( zji,zjj,jl+1) .GT. epsi10 ) .AND. &317 ( ht_i( zji,zjj,jl+1).LE. zhbnew(zji,zjj,jl) ) &352 zremap_flag(ii,ij) = 0 353 ELSEIF ( ( a_i(ii,ij,jl+1) .GT. epsi10 ) .AND. & 354 ( ht_i(ii,ij,jl+1).LE. zhbnew(ii,ij,jl) ) & 318 355 ) THEN 319 zremap_flag( zji,zjj) = 0356 zremap_flag(ii,ij) = 0 320 357 ENDIF 321 358 322 359 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 323 360 ! jl, ji 324 IF (zhbnew( zji,zjj,jl).gt.hi_max(jl+1)) THEN325 zremap_flag( zji,zjj) = 0361 IF (zhbnew(ii,ij,jl).gt.hi_max(jl+1)) THEN 362 zremap_flag(ii,ij) = 0 326 363 ENDIF 327 364 ! jl, ji 328 IF (zhbnew( zji,zjj,jl).lt.hi_max(jl-1)) THEN329 zremap_flag( zji,zjj) = 0365 IF (zhbnew(ii,ij,jl).lt.hi_max(jl-1)) THEN 366 zremap_flag(ii,ij) = 0 330 367 ENDIF 331 368 ! jl, ji … … 379 416 !- 7.2 Area lost due to melting of thin ice (first category, klbnd) 380 417 DO ji = 1, nbrem 381 zji = nind_i(ji)382 zjj = nind_j(ji)418 ii = nind_i(ji) 419 ij = nind_j(ji) 383 420 384 421 !ji 385 IF (a_i( zji,zjj,klbnd) .gt. epsi10) THEN386 zdh0 = zdhice( zji,zjj,klbnd) !decrease of ice thickness in the lower category422 IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN 423 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 387 424 ! ji, a_i > epsi10 388 425 IF (zdh0 .lt. 0.0) THEN !remove area from category 1 … … 391 428 392 429 !Integrate g(1) from 0 to dh0 to estimate area melted 393 zetamax = MIN(zdh0,hR( zji,zjj,klbnd)) - hL(zji,zjj,klbnd)430 zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd) 394 431 IF (zetamax.gt.0.0) THEN 395 432 zx1 = zetamax 396 433 zx2 = 0.5 * zetamax*zetamax 397 zda0 = g1( zji,zjj,klbnd) * zx2 + g0(zji,zjj,klbnd) * zx1 !ice area removed434 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 398 435 ! Constrain new thickness <= ht_i 399 zdamax = a_i( zji,zjj,klbnd) * &400 (1.0 - ht_i( zji,zjj,klbnd)/zht_i_o(zji,zjj,klbnd)) ! zdamax > 0436 zdamax = a_i(ii,ij,klbnd) * & 437 (1.0 - ht_i(ii,ij,klbnd)/zht_i_o(ii,ij,klbnd)) ! zdamax > 0 401 438 !ice area lost due to melting of thin ice 402 439 zda0 = MIN(zda0, zdamax) 403 440 404 441 ! Remove area, conserving volume 405 ht_i( zji,zjj,klbnd) = ht_i(zji,zjj,klbnd) &406 * a_i( zji,zjj,klbnd) / ( a_i(zji,zjj,klbnd) - zda0 )407 a_i( zji,zjj,klbnd) = a_i(zji,zjj,klbnd) - zda0408 v_i( zji,zjj,klbnd) = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd)442 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 443 * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 444 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 445 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem@useless ? 409 446 ENDIF ! zetamax > 0 410 447 ! ji, a_i > epsi10 … … 412 449 ELSE ! if ice accretion 413 450 ! ji, a_i > epsi10; zdh0 > 0 414 IF ( ntyp .EQ. 1 ) zhbnew( zji,zjj,klbnd-1) = MIN(zdh0,hi_max(klbnd))451 IF ( ntyp .EQ. 1 ) zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 415 452 ! zhbnew was 0, and is shifted to the right to account for thin ice 416 453 ! growth in openwater (F0 = f1) 417 IF ( ntyp .NE. 1 ) zhbnew( zji,zjj,0) = 0454 IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0 418 455 ! in other types there is 419 456 ! no open water growth (F0 = 0) … … 446 483 447 484 DO ji = 1, nbrem 448 zji = nind_i(ji)449 zjj = nind_j(ji)450 451 IF (zhbnew( zji,zjj,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1485 ii = nind_i(ji) 486 ij = nind_j(ji) 487 488 IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 452 489 453 490 ! left and right integration limits in eta space 454 zvetamin(ji) = MAX(hi_max(jl), hL( zji,zjj,jl)) - hL(zji,zjj,jl)455 zvetamax(ji) = MIN(zhbnew( zji,zjj,jl), hR(zji,zjj,jl)) - hL(zji,zjj,jl)456 zdonor( zji,zjj,jl) = jl491 zvetamin(ji) = MAX(hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl) 492 zvetamax(ji) = MIN(zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl) 493 zdonor(ii,ij,jl) = jl 457 494 458 495 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl … … 460 497 ! left and right integration limits in eta space 461 498 zvetamin(ji) = 0.0 462 zvetamax(ji) = MIN(hi_max(jl), hR( zji,zjj,jl+1)) - hL(zji,zjj,jl+1)463 zdonor( zji,zjj,jl) = jl + 1499 zvetamax(ji) = MIN(hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1) 500 zdonor(ii,ij,jl) = jl + 1 464 501 465 502 ENDIF ! zhbnew(jl) > hi_max(jl) … … 475 512 zwk2 = zwk2 * zetamax 476 513 zx3 = 1.0/3.0 * (zwk2 - zwk1) 477 nd = zdonor( zji,zjj,jl)478 zdaice( zji,zjj,jl) = g1(zji,zjj,nd)*zx2 + g0(zji,zjj,nd)*zx1479 zdvice( zji,zjj,jl) = g1(zji,zjj,nd)*zx3 + g0(zji,zjj,nd)*zx2 + &480 zdaice( zji,zjj,jl)*hL(zji,zjj,nd)514 nd = zdonor(ii,ij,jl) 515 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 516 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + & 517 zdaice(ii,ij,jl)*hL(ii,ij,nd) 481 518 482 519 END DO ! ji … … 493 530 494 531 DO ji = 1, nbrem 495 zji = nind_i(ji) 496 zjj = nind_j(ji) 497 IF ( ( zhimin .GT. 0.0 ) .AND. & 498 ( ( a_i(zji,zjj,1) .GT. epsi10 ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) & 499 ) THEN 500 a_i(zji,zjj,1) = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin 501 ht_i(zji,zjj,1) = zhimin 502 v_i(zji,zjj,1) = a_i(zji,zjj,1)*ht_i(zji,zjj,1) 532 ii = nind_i(ji) 533 ij = nind_j(ji) 534 IF ( ( a_i(ii,ij,1) > epsi10 ) .AND. ( ht_i(ii,ij,1) < hiclim ) ) THEN 535 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim 536 ht_i(ii,ij,1) = hiclim 537 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem@useless 503 538 ENDIF 504 539 END DO !ji … … 625 660 626 661 INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices 627 INTEGER :: zji, zjj ! indices when changing from 2D-1D is done662 INTEGER :: ii, ij ! indices when changing from 2D-1D is done 628 663 629 664 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn … … 759 794 760 795 DO ji = 1, nbrem 761 zji = nind_i(ji)762 zjj = nind_j(ji)763 764 jl1 = zdonor( zji,zjj,jl)765 zindb = MAX( 0.0 , SIGN( 1.0 , v_i( zji,zjj,jl1) - epsi10 ) )766 zworka( zji,zjj) = zdvice(zji,zjj,jl) / MAX(v_i(zji,zjj,jl1),epsi10) * zindb796 ii = nind_i(ji) 797 ij = nind_j(ji) 798 799 jl1 = zdonor(ii,ij,jl) 800 zindb = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 801 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * zindb 767 802 IF( jl1 == jl) THEN ; jl2 = jl1+1 768 803 ELSE ; jl2 = jl … … 773 808 !-------------- 774 809 775 a_i( zji,zjj,jl1) = a_i(zji,zjj,jl1) - zdaice(zji,zjj,jl)776 a_i( zji,zjj,jl2) = a_i(zji,zjj,jl2) + zdaice(zji,zjj,jl)810 a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 811 a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 777 812 778 813 !-------------- … … 780 815 !-------------- 781 816 782 v_i( zji,zjj,jl1) = v_i(zji,zjj,jl1) - zdvice(zji,zjj,jl)783 v_i( zji,zjj,jl2) = v_i(zji,zjj,jl2) + zdvice(zji,zjj,jl)817 v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 818 v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 784 819 785 820 !-------------- … … 787 822 !-------------- 788 823 789 zdvsnow = v_s( zji,zjj,jl1) * zworka(zji,zjj)790 v_s( zji,zjj,jl1) = v_s(zji,zjj,jl1) - zdvsnow791 v_s( zji,zjj,jl2) = v_s(zji,zjj,jl2) + zdvsnow824 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 825 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 826 v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow 792 827 793 828 !-------------------- … … 795 830 !-------------------- 796 831 797 zdesnow = e_s( zji,zjj,1,jl1) * zworka(zji,zjj)798 e_s( zji,zjj,1,jl1) = e_s(zji,zjj,1,jl1) - zdesnow799 e_s( zji,zjj,1,jl2) = e_s(zji,zjj,1,jl2) + zdesnow832 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 833 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow 834 e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow 800 835 801 836 !-------------- … … 803 838 !-------------- 804 839 805 zdo_aice = oa_i( zji,zjj,jl1) * zdaice(zji,zjj,jl)806 oa_i( zji,zjj,jl1) = oa_i(zji,zjj,jl1) - zdo_aice807 oa_i( zji,zjj,jl2) = oa_i(zji,zjj,jl2) + zdo_aice840 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 841 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice 842 oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice 808 843 809 844 !-------------- … … 811 846 !-------------- 812 847 813 zdsm_vice = smv_i( zji,zjj,jl1) * zworka(zji,zjj)814 smv_i( zji,zjj,jl1) = smv_i(zji,zjj,jl1) - zdsm_vice815 smv_i( zji,zjj,jl2) = smv_i(zji,zjj,jl2) + zdsm_vice848 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 849 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice 850 smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice 816 851 817 852 !--------------------- … … 819 854 !--------------------- 820 855 821 zdaTsf = t_su( zji,zjj,jl1) * zdaice(zji,zjj,jl)822 zaTsfn( zji,zjj,jl1) = zaTsfn(zji,zjj,jl1) - zdaTsf823 zaTsfn( zji,zjj,jl2) = zaTsfn(zji,zjj,jl2) + zdaTsf856 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 857 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf 858 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 824 859 825 860 END DO ! ji … … 832 867 !CDIR NODEP 833 868 DO ji = 1, nbrem 834 zji = nind_i(ji)835 zjj = nind_j(ji)836 837 jl1 = zdonor( zji,zjj,jl)869 ii = nind_i(ji) 870 ij = nind_j(ji) 871 872 jl1 = zdonor(ii,ij,jl) 838 873 IF (jl1 .EQ. jl) THEN 839 874 jl2 = jl+1 … … 842 877 ENDIF 843 878 844 zdeice = e_i( zji,zjj,jk,jl1) * zworka(zji,zjj)845 e_i( zji,zjj,jk,jl1) = e_i(zji,zjj,jk,jl1) - zdeice846 e_i( zji,zjj,jk,jl2) = e_i(zji,zjj,jk,jl2) + zdeice879 zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij) 880 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice 881 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 847 882 END DO ! ji 848 883 END DO ! jk … … 860 895 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 861 896 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 862 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl) )) !0 if no ice and 1 if yes897 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 863 898 ELSE 864 899 ht_i(ji,jj,jl) = 0._wp … … 967 1002 zshiftflag = 1 968 1003 zdonor(ji,jj,jl) = jl 969 zdaice(ji,jj,jl) = a_i(ji,jj,jl) 970 zdvice(ji,jj,jl) = v_i(ji,jj,jl) 1004 ! begin TECLIM change 1005 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) 1006 !zdvice(ji,jj,jl) = v_i(ji,jj,jl) 1007 zdaice(ji,jj,jl) = a_i(ji,jj,jl)/2 1008 zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1))/2 1009 ! end TECLIM change 971 1010 ENDIF 972 1011 END DO ! ji -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90
r3625 r4045 26 26 27 27 !!---------------------------------------------------------------------- 28 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)28 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 29 29 !! $Id$ 30 30 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r3791 r4045 41 41 USE agrif_lim2_interp 42 42 #endif 43 #if defined key_bdy 44 USE bdyice_lim 45 #endif 43 46 44 47 IMPLICIT NONE … … 53 56 # include "vectopt_loop_substitute.h90" 54 57 !!---------------------------------------------------------------------- 55 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)58 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 56 59 !! $Id$ 57 60 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 473 476 CALL lbc_lnk( zs12(:,:), 'F', 1. ) 474 477 478 !#if defined key_bdy 479 ! ! clem: change zs1, zs2, zs12 at the boundary for each iteration 480 ! CALL bdy_ice_lim_dyn( 2, zs1, zs2, zs12 ) 481 ! CALL lbc_lnk( zs1 (:,:), 'T', 1. ) 482 ! CALL lbc_lnk( zs2 (:,:), 'T', 1. ) 483 ! CALL lbc_lnk( zs12(:,:), 'F', 1. ) 484 !#endif 485 475 486 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 476 487 DO jj = k_j1+1, k_jpj-1 … … 520 531 521 532 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 522 #if defined key_agrif 533 #if defined key_agrif && defined key_lim2 523 534 CALL agrif_rhg_lim2( jter, nevp, 'U' ) 524 535 #endif … … 548 559 549 560 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 550 #if defined key_agrif 561 #if defined key_agrif && defined key_lim2 551 562 CALL agrif_rhg_lim2( jter, nevp, 'V' ) 552 563 #endif … … 577 588 578 589 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 579 #if defined key_agrif 590 #if defined key_agrif && defined key_lim2 580 591 CALL agrif_rhg_lim2( jter, nevp , 'V' ) 581 592 #endif … … 608 619 609 620 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 610 #if defined key_agrif 621 #if defined key_agrif && defined key_lim2 611 622 CALL agrif_rhg_lim2( jter, nevp, 'U' ) 612 623 #endif 613 624 614 625 ENDIF 626 627 !#if defined key_bdy 628 ! ! clem: change u_ice and v_ice at the boundary for each iteration 629 ! CALL bdy_ice_lim_dyn( 1 ) 630 !#endif 615 631 616 632 IF(ln_ctl) THEN … … 624 640 ENDIF 625 641 626 ! 642 ! ! ==================== ! 627 643 END DO ! end loop over jter ! 628 644 ! ! ==================== ! 629 630 645 ! 631 646 !------------------------------------------------------------------------------! 632 647 ! 4) Prevent ice velocities when the ice is thin 633 648 !------------------------------------------------------------------------------! 634 ! 635 ! If the ice thickness is below 1cm then ice velocity should equal the 649 !clem : add hminrhg in the namelist 650 ! 651 ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 636 652 ! ocean velocity, 637 653 ! This prevents high velocity when ice is thin … … 642 658 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 643 659 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 644 IF ( zdummy .LE. 5.0e-2 ) THEN 660 !zdummy = vt_i(ji,jj) 661 IF ( zdummy .LE. hminrhg ) THEN 645 662 u_ice(ji,jj) = u_oce(ji,jj) 646 663 v_ice(ji,jj) = v_oce(ji,jj) … … 651 668 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 652 669 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 653 #if defined key_agrif 670 #if defined key_agrif && defined key_lim2 654 671 CALL agrif_rhg_lim2( nevp , nevp, 'U' ) 655 672 CALL agrif_rhg_lim2( nevp , nevp, 'V' ) 656 673 #endif 674 #if defined key_bdy 675 ! clem: change u_ice and v_ice at the boundary 676 CALL bdy_ice_lim_dyn( 1 ) 677 #endif 657 678 658 679 DO jj = k_j1+1, k_jpj-1 … … 660 681 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 661 682 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 662 IF ( zdummy .LE. 5.0e-2 ) THEN 683 !zdummy = vt_i(ji,jj) 684 IF ( zdummy .LE. hminrhg ) THEN 663 685 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 664 686 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & … … 684 706 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 685 707 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 686 687 IF ( zdummy .LE. 5.0e-2) THEN708 !zdummy = vt_i(ji,jj) 709 IF ( zdummy .LE. hminrhg ) THEN 688 710 689 711 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & … … 738 760 divu_i (ji,jj) = zdd (ji,jj) 739 761 delta_i(ji,jj) = deltat(ji,jj) 762 ! begin TECLIM change 763 zdst(ji,jj)= ( e2u( ji , jj ) * v_ice1(ji,jj) & 764 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 765 & + e1v( ji , jj ) * u_ice2(ji,jj) & 766 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 740 767 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst(ji,jj) * zdst(ji,jj) ) 768 ! end TECLIM change 741 769 END DO 742 770 END DO 743 CALL lbc_lnk( divu_i (:,:), 'T', 1. ) ! Lateral boundary condition 771 772 ! Lateral boundary condition 773 CALL lbc_lnk( divu_i (:,:), 'T', 1. ) 744 774 CALL lbc_lnk( delta_i(:,:), 'T', 1. ) 775 ! CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 745 776 CALL lbc_lnk( shear_i(:,:), 'T', 1. ) 746 777 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
r3625 r4045 38 38 39 39 !!---------------------------------------------------------------------- 40 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)40 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 41 41 !! $Id$ 42 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 162 162 CALL iom_rstput( iter, nitrst, numriw, 'v_ice' , v_ice ) 163 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 164 166 CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) 165 167 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) … … 340 342 !Control of date 341 343 342 IF( ( nit000 - INT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 ) &344 IF( ( nit000 - NINT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 ) & 343 345 & CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nit000 in ice restart', & 344 346 & ' verify the file or rerun with the value 0 for the', & 345 347 & ' control of time parameter nrstdt' ) 346 IF( INT(zfice) /= nn_fsbc .AND. ABS( nrstdt ) == 1 ) &348 IF( NINT(zfice) /= nn_fsbc .AND. ABS( nrstdt ) == 1 ) & 347 349 & CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nn_fsbc in ice restart', & 348 350 & ' verify the file or rerun with the value 0 for the', & … … 437 439 CALL iom_get( numrir, jpdom_autoglo, 'v_ice' , v_ice ) 438 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 439 443 CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i ) 440 444 CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i ) … … 563 567 END DO 564 568 ! 565 CALL iom_close( numrir )569 !clem CALL iom_close( numrir ) 566 570 ! 567 571 CALL wrk_dealloc( nlay_i, zs_zero ) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r3905 r4045 10 10 !! ! + simplification of the ice-ocean stress calculation 11 11 !! 3.4 ! 2011-02 (G. Madec) dynamical allocation 12 !! - ! 2012 (D. Iovino) salt flux change 13 !! - ! 2012-05 (C. Rousset) add penetration solar flux 12 14 !! 3.5 ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass 13 15 !!---------------------------------------------------------------------- … … 35 37 USE prtctl ! Print control 36 38 USE cpl_oasis3, ONLY : lk_cpl 39 USE traqsr ! clem: add penetration of solar flux into the calculation of heat budget 37 40 USE oce, ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n 38 41 USE dom_ice, ONLY : tms … … 57 60 # include "vectopt_loop_substitute.h90" 58 61 !!---------------------------------------------------------------------- 59 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)62 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 60 63 !! $Id$ 61 64 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 106 109 REAL(wp) :: zfcm1 , zfcm2 ! - - 107 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb, zalbp ! 2D/3D workspace 111 REAL(wp) :: zzfcm1, zfscmbq ! clem: for light penetration 108 112 !!--------------------------------------------------------------------- 109 113 … … 119 123 DO ji = 1, jpi 120 124 zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 121 ifvt = zinda * MAX( rzero , SIGN( rone, - phicif(ji,jj) ) ) !subscripts are bad here122 i1mfr = 1.0 - MAX( rzero , SIGN( rone , - ( at_i(ji,jj)) ) )125 ifvt = zinda * MAX( rzero , SIGN( rone, - phicif(ji,jj) ) ) !subscripts are bad here 126 i1mfr = 1.0 - MAX( rzero , SIGN( rone , - at_i(ji,jj) ) ) 123 127 idfr = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) ) 124 128 iflt = zinda * (1 - i1mfr) * (1 - ifvt ) … … 141 145 142 146 ! computation the solar flux at ocean surface 143 zfcm1 = pfrld(ji,jj) * qsr(ji,jj) + ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) 147 zfcm1 = pfrld(ji,jj) * qsr(ji,jj) + & 148 & ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 144 149 ! fstric Solar flux transmitted trough the ice 145 150 ! qsr Net short wave heat flux on free ocean 146 151 ! new line 147 fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) 152 fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 153 154 ! solar flux and fscmbq with light penetration (clem) 155 zzfcm1 = pfrld(ji,jj) * qsr(ji,jj) * oatte(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) 156 zfscmbq = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) 148 157 149 158 ! computation the non solar heat flux at ocean surface 150 zfcm2 = - z fcm1 & ! ???151 & + iflt * fscmbq(ji,jj)& ! total ablation: heat given to the ocean159 zfcm2 = - zzfcm1 & ! 160 & + iflt * zfscmbq & ! total ablation: heat given to the ocean 152 161 & + ifral * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice & 153 162 & + ifrdv * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * r1_rdtice & … … 170 179 ! ! fdtcn : turbulent oceanic heat flux 171 180 172 !!gm this IF prevents the vertorisation of the whole loop173 IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN174 WRITE(numout,*) ' lim_sbc : heat fluxes '175 WRITE(numout,*) ' qsr : ', qsr(jiindx,jjindx)176 WRITE(numout,*) ' pfrld : ', pfrld(jiindx,jjindx)177 WRITE(numout,*) ' fstric : ', fstric (jiindx,jjindx)178 WRITE(numout,*)179 WRITE(numout,*) ' qns : ', qns(jiindx,jjindx)180 WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)181 WRITE(numout,*) ' ifral : ', ifral182 WRITE(numout,*) ' ial : ', ial183 WRITE(numout,*) ' qcmif : ', qcmif(jiindx,jjindx)184 WRITE(numout,*) ' qldif : ', qldif(jiindx,jjindx)185 186 187 WRITE(numout,*) ' ifrdv : ', ifrdv188 WRITE(numout,*) ' qfvbq : ', qfvbq(jiindx,jjindx)189 WRITE(numout,*) ' qdtcn : ', qdtcn(jiindx,jjindx)190 191 192 WRITE(numout,*) ' '193 WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)194 WRITE(numout,*) ' fhmec : ', fhmec(jiindx,jjindx)195 WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx)196 WRITE(numout,*) ' fhbri : ', fhbri(jiindx,jjindx)197 WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)198 ENDIF199 !!gm end181 !!gm this IF prevents the vertorisation of the whole loop 182 ! IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 183 ! WRITE(numout,*) ' lim_sbc : heat fluxes ' 184 ! WRITE(numout,*) ' qsr : ', qsr(jiindx,jjindx) 185 ! WRITE(numout,*) ' pfrld : ', pfrld(jiindx,jjindx) 186 ! WRITE(numout,*) ' fstric : ', fstric (jiindx,jjindx) 187 ! WRITE(numout,*) 188 ! WRITE(numout,*) ' qns : ', qns(jiindx,jjindx) 189 ! WRITE(numout,*) ' fdtcn : ', fdtcn(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_mec : ', fheat_mec(jiindx,jjindx) 205 ! WRITE(numout,*) ' fhbri : ', fhbri(jiindx,jjindx) 206 ! WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 207 ! ENDIF 208 !!gm end 200 209 END DO 201 210 END DO … … 370 379 !! ** input : Namelist namicedia 371 380 !!------------------------------------------------------------------- 381 REAL(wp) :: zsum, zarea 372 382 ! 373 383 INTEGER :: ji, jj ! dummy loop indices … … 390 400 END WHERE 391 401 ENDIF 402 ! clem modif 403 iatte(:,:) = 1._wp 404 oatte(:,:) = 1._wp 405 ! 392 406 ! ! embedded sea ice 393 407 IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass … … 435 449 ENDIF 436 450 ! 451 !!? IF( .NOT. ln_rstart ) THEN ! delete the initial ssh below sea-ice area 452 !!? ! 453 !!? zarea = glob_sum( e1e2t(:,:) ) ! interior global domain surface 454 !!? zsum = glob_sum( e1e2t(:,:) * ( snwice_mass(:,:) ) ) / zarea * r1_rau0 455 !!? sshn(:,:) = sshn(:,:) - zsum 456 !!? sshb(:,:) = sshb(:,:) - zsum 457 !!? ENDIF 458 ! 459 437 460 END SUBROUTINE lim_sbc_init 438 461 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90
r3625 r4045 20 20 21 21 !!---------------------------------------------------------------------- 22 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2010)22 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 23 23 !! $Id$ 24 24 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r3625 r4045 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 … … 92 93 REAL(wp) :: zfntlat, zpareff, zareamin, zcoef ! - - 93 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) 96 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 94 97 !!------------------------------------------------------------------- 95 98 96 99 CALL wrk_alloc( jpi, jpj, zqlbsbq ) 97 100 101 ! ------------------------------- 102 !- check conservation (C Rousset) 103 IF (ln_limdiahsb) THEN 104 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 105 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 106 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 107 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 108 ENDIF 109 !- check conservation (C Rousset) 110 ! ------------------------------- 111 98 112 !------------------------------------------------------------------------------! 99 113 ! 1) Initialization of diagnostic variables ! … … 109 123 DO ji = 1, jpi 110 124 !Energy of melting q(S,T) [J.m-3] 111 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * nlay_i125 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 ) 112 126 !0 if no ice and 1 if yes 113 127 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) ) ) … … 121 135 DO ji = 1, jpi 122 136 !Energy of melting q(S,T) [J.m-3] 123 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * nlay_s137 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 ) 124 138 !0 if no ice and 1 if yes 125 139 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) ) ) … … 134 148 ! 1.3) Set some dummies to 0 135 149 !----------------------------- 136 rdvosif(:,:) = 0.e0 ! variation of ice volume at surface137 rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom138 fdvolif(:,:) = 0.e0 ! total variation of ice volume139 rdvonif(:,:) = 0.e0 ! lateral variation of ice volume140 fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice141 ffltbif(:,:) = 0.e0 ! linked with fstric142 qfvbq (:,:) = 0.e0 ! linked with fstric143 rdm_snw(:,:) = 0.e0 ! variation of snow mass per unit area144 rdm_ice(:,:) = 0.e0 ! variation of ice mass per unit area145 hicifp (:,:) = 0.e0 ! daily thermodynamic ice production.146 sfx_bri(:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean147 fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean148 sfx_thd(:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay150 !clem rdvosif(:,:) = 0.e0 ! variation of ice volume at surface 151 !clem rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom 152 !clem fdvolif(:,:) = 0.e0 ! total variation of ice volume 153 !clem rdvonif(:,:) = 0.e0 ! lateral variation of ice volume 154 !clem fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice 155 !clem ffltbif(:,:) = 0.e0 ! linked with fstric 156 !clem qfvbq (:,:) = 0.e0 ! linked with fstric 157 !clem rdm_snw(:,:) = 0.e0 ! variation of snow mass per unit area 158 !clem rdm_ice(:,:) = 0.e0 ! variation of ice mass per unit area 159 !clem hicifp (:,:) = 0.e0 ! daily thermodynamic ice production. 160 !clem sfx_bri(:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean 161 !clem fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean 162 !clem sfx_thd(:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay 149 163 150 164 !----------------------------------- … … 165 179 !CDIR NOVERRCHK 166 180 DO ji = 1, jpi 167 zthsnice = SUM( ht_s(ji,jj,1:jpl) ) + SUM( ht_i(ji,jj,1:jpl) )168 zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) )169 181 phicif(ji,jj) = vt_i(ji,jj) 170 182 pfrld(ji,jj) = 1.0 - at_i(ji,jj) 171 zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )183 zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) ) 172 184 ! 173 185 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 180 192 181 193 ! here the drag will depend on ice thickness and type (0.006) 182 fdtcn(ji,jj) = zind b* rau0 * rcp * 0.006 * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) )194 fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) ) 183 195 ! also category dependent 184 196 ! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead 185 qdtcn(ji,jj) = zind b* fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice197 qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 186 198 ! 187 199 ! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) 188 200 ! ! caution: exponent betas used as more snow can fallinto leads 189 201 qldif(ji,jj) = tms(ji,jj) * rdt_ice * ( & 190 & pfrld(ji,jj) * ( qsr(ji,jj) & ! solar heat202 & pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif 191 203 & + qns(ji,jj) & ! non solar heat 192 204 & + fdtcn(ji,jj) & ! turbulent ice-ocean heat 193 & + fsbbq(ji,jj) * ( 1.0 - zind b) ) & ! residual heat from previous step205 & + fsbbq(ji,jj) * ( 1.0 - zinda ) ) & ! residual heat from previous step 194 206 & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus ) ! latent heat of sprecip melting 195 207 ! … … 206 218 ! 207 219 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 208 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda )220 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) 209 221 ! 210 222 ! oceanic heat flux (limthd_dh) 211 fbif (ji,jj) = zind b* ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) )223 fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 212 224 ! 213 225 END DO … … 294 306 CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb), qfvbq , jpi, jpj, npb(1:nbpb) ) 295 307 308 CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) ) ! clem modif 309 CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) ) ! clem modif 296 310 !-------------------------------- 297 311 ! 4.3) Thermodynamic processes … … 411 425 ! 5.4) Diagnostic thermodynamic growth rates 412 426 !-------------------------------------------- 413 d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes414 427 !clem@useless d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes 428 !clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 415 429 416 430 IF( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) … … 448 462 ENDIF 449 463 ! 464 ! ------------------------------- 465 !- check conservation (C Rousset) 466 IF (ln_limdiahsb) THEN 467 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 468 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 469 470 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 471 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 472 473 zchk_vmin = glob_min(v_i) 474 zchk_amax = glob_max(SUM(a_i,dim=3)) 475 zchk_amin = glob_min(a_i) 476 477 IF(lwp) THEN 478 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limthd) = ',(zchk_v_i * rday) 479 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * rday) 480 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limthd) = ',(zchk_vmin * 1.e-3) 481 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limthd) = ',zchk_amax 482 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limthd) = ',zchk_amin 483 ENDIF 484 ENDIF 485 !- check conservation (C Rousset) 486 ! ------------------------------- 487 ! 450 488 CALL wrk_dealloc( jpi, jpj, zqlbsbq ) 451 489 ! … … 472 510 DO jk = 1, nlay_i ! total q over all layers, ice [J.m-2] 473 511 DO ji = kideb, kiut 474 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 ) 475 513 eti (ji,jl) = eti(ji,jl) + etilayer(ji,jk) 476 514 END DO 477 515 END DO 478 516 DO ji = kideb, kiut ! total q over all layers, snow [J.m-2] 479 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 ) 480 518 END DO 481 519 ! … … 498 536 499 537 INTEGER :: ji, jk ! loop indices 500 INTEGER :: zji, zjj538 INTEGER :: ii, ij 501 539 INTEGER :: numce ! number of points for which conservation is violated 502 540 REAL(wp) :: meance ! mean conservation error … … 521 559 !---------------------------------------- 522 560 DO ji = kideb, kiut 523 zji = MOD( npb(ji) - 1 , jpi ) + 1524 zjj = ( npb(ji) - 1 ) / jpi + 1561 ii = MOD( npb(ji) - 1 , jpi ) + 1 562 ij = ( npb(ji) - 1 ) / jpi + 1 525 563 fatm (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) 526 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc( zji,zjj,jl)564 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(ii,ij,jl) 527 565 END DO 528 566 … … 579 617 IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. & 580 618 ( cons_error(ji,jl) .GT. max_cons_err ) ) THEN 581 zji = MOD( npb(ji) - 1, jpi ) + 1582 zjj = ( npb(ji) - 1 ) / jpi + 1619 ii = MOD( npb(ji) - 1, jpi ) + 1 620 ij = ( npb(ji) - 1 ) / jpi + 1 583 621 ! 584 622 WRITE(numout,*) ' alerte 1 ' … … 586 624 WRITE(numout,*) ' heat diffusion in the ice ' 587 625 WRITE(numout,*) ' Category : ', jl 588 WRITE(numout,*) ' zji , zjj : ', zji, zjj589 WRITE(numout,*) ' lat, lon : ', gphit( zji,zjj), glamt(zji,zjj)626 WRITE(numout,*) ' ii , ij : ', ii, ij 627 WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij) 590 628 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 591 629 WRITE(numout,*) ' surf_error : ', surf_error(ji,jl) … … 615 653 WRITE(numout,*) ' fc_bo : ', - fc_bo_i (ji) 616 654 WRITE(numout,*) ' foc : ', fbif_1d(ji) 617 WRITE(numout,*) ' fstroc : ', fstroc ( zji,zjj,jl)655 WRITE(numout,*) ' fstroc : ', fstroc (ii,ij,jl) 618 656 WRITE(numout,*) ' i0 : ', i0(ji) 619 657 WRITE(numout,*) ' qsr_ice : ', (1.0-i0(ji))*qsr_ice_1d(ji) … … 651 689 ! 652 690 INTEGER :: ji ! loop indices 653 INTEGER :: zji, zjj, numce ! local integers691 INTEGER :: ii, ij, numce ! local integers 654 692 REAL(wp) :: meance, max_cons_err !local scalar 655 693 !!--------------------------------------------------------------------- … … 669 707 !---------------------------------------- 670 708 DO ji = kideb, kiut 671 zji = MOD( npb(ji) - 1 , jpi ) + 1672 zjj = ( npb(ji) - 1 ) / jpi + 1709 ii = MOD( npb(ji) - 1 , jpi ) + 1 710 ij = ( npb(ji) - 1 ) / jpi + 1 673 711 674 712 fatm (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji) ! total heat flux 675 sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc( zji,zjj,jl)713 sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(ii,ij,jl) 676 714 cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 677 715 END DO … … 706 744 DO ji = kideb, kiut 707 745 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 708 zji = MOD( npb(ji) - 1, jpi ) + 1709 zjj = ( npb(ji) - 1 ) / jpi + 1746 ii = MOD( npb(ji) - 1, jpi ) + 1 747 ij = ( npb(ji) - 1 ) / jpi + 1 710 748 ! 711 749 WRITE(numout,*) ' alerte 1 - category : ', jl 712 750 WRITE(numout,*) ' Untolerated conservation error after limthd_ent ' 713 WRITE(numout,*) ' zji , zjj : ', zji, zjj714 WRITE(numout,*) ' lat, lon : ', gphit( zji,zjj), glamt(zji,zjj)751 WRITE(numout,*) ' ii , ij : ', ii, ij 752 WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij) 715 753 WRITE(numout,*) ' * ' 716 754 WRITE(numout,*) ' Ftotal : ', sum_fluxq(ji,jl) … … 724 762 WRITE(numout,*) ' foce : ', fbif_1d(ji) 725 763 WRITE(numout,*) ' fres : ', ftotal_fin(ji) 726 WRITE(numout,*) ' fhbri : ', fhbricat( zji,zjj,jl)764 WRITE(numout,*) ' fhbri : ', fhbricat(ii,ij,jl) 727 765 WRITE(numout,*) ' * ' 728 766 WRITE(numout,*) ' Heat contents --- : ' … … 792 830 !!------------------------------------------------------------------- 793 831 NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, & 794 & hicmin, hiclim, amax ,&832 & hicmin, hiclim, & 795 833 & sbeta , parlat, hakspl, hibspl, exld, & 796 834 & hakdif, hnzst , thth , parsub, alphs, betas, & … … 818 856 WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin 819 857 WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim 820 WRITE(numout,*)' maximum lead fraction amax = ', amax821 858 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 822 859 WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r3808 r4045 6 6 !! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D 7 7 !! ! 2005-06 (M. Vancoppenolle) 3D version 8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm snif & rdmicif8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw & rdm_ice 9 9 !! 3.4 ! 2011-02 (G. Madec) dynamical allocation 10 10 !! 3.5 ! 2012-10 (G. Madec & co) salt flux + bug fixes … … 39 39 40 40 !!---------------------------------------------------------------------- 41 !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)41 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 42 42 !! $Id$ 43 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 73 73 !! 74 74 INTEGER :: ji , jk ! dummy loop indices 75 INTEGER :: ii, ij ! 2D corresponding indices to ji75 INTEGER :: ii, ij ! 2D corresponding indices to ji 76 76 INTEGER :: isnow ! switch for presence (1) or absence (0) of snow 77 77 INTEGER :: isnowic ! snow ice formation not … … 107 107 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_pre ! snow precipitation 108 108 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_sub ! snow sublimation 109 REAL(wp), POINTER, DIMENSION(:) :: zsfx_melt ! salt flux due to ice melt110 109 111 110 REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah … … 120 119 REAL(wp), POINTER, DIMENSION(:,:) :: zqt_i_lay ! total ice heat content 121 120 121 ! mass and salt flux (clem) 122 REAL(wp) :: zdvres, zdvsur, zdvbot 123 REAL(wp), POINTER, DIMENSION(:) :: zviold, zvsold ! old ice volume... 124 122 125 ! Heat conservation 123 126 INTEGER :: num_iter_max, numce_dh … … 128 131 129 132 CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 130 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, z sfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )133 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 131 134 CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i ) 132 135 CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) 133 136 134 zsfx_melt (:) = 0._wp 137 CALL wrk_alloc( jpij, zviold, zvsold ) ! clem 138 135 139 ftotal_fin(:) = 0._wp 136 140 zfdt_init (:) = 0._wp 137 141 zfdt_final(:) = 0._wp 138 142 143 dh_i_surf (:) = 0._wp 144 dh_i_bott (:) = 0._wp 145 dh_snowice(:) = 0._wp 146 139 147 DO ji = kideb, kiut 140 148 old_ht_i_b(ji) = ht_i_b(ji) 141 149 old_ht_s_b(ji) = ht_s_b(ji) 150 zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem 151 zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem 142 152 END DO 143 153 ! … … 164 174 ! 165 175 DO ji = kideb, kiut ! Layer thickness 166 zh_i(ji) = ht_i_b(ji) / nlay_i167 zh_s(ji) = ht_s_b(ji) / nlay_s176 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 177 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 168 178 END DO 169 179 ! … … 171 181 DO jk = 1, nlay_s 172 182 DO ji = kideb, kiut 173 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s183 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) 174 184 END DO 175 185 END DO … … 178 188 DO jk = 1, nlay_i 179 189 DO ji = kideb, kiut 180 zzc = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i190 zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 181 191 zqt_i(ji) = zqt_i(ji) + zzc 182 192 zqt_i_lay(ji,jk) = zzc … … 244 254 zhn = 1.0 - MAX( zzero , SIGN( zone , - zhsnew ) ) 245 255 ht_s_b(ji) = MAX( zzero , zhsnew ) 256 ! we recompute dh_s_tot (clem) 257 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 246 258 ! Volume and mass variations of snow 247 259 dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) 248 260 dvsbq_1d (ji) = MIN( zzero, dvsbq_1d(ji) ) 249 rdm_snw_1d(ji) = rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji)261 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 250 262 END DO ! ji 251 263 … … 254 266 !-------------------------- 255 267 DO ji = kideb, kiut 256 dh_i_surf(ji) = 0._wp257 268 z_f_surf (ji) = zqfont_su(ji) * r1_rdtice ! heat conservation test 258 269 zdq_i (ji) = 0._wp … … 272 283 zdq_i (ji ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 273 284 ! 274 ! ! contribution to ice-ocean salt flux 275 zsfx_melt(ji) = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 285 ! clem 286 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 287 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 276 288 END DO 277 289 END DO … … 334 346 DO ji = kideb,kiut 335 347 q_s_b (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 336 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s! heat conservation348 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) ! heat conservation 337 349 END DO 338 350 END DO … … 375 387 ! Basal growth rate = - F*dt / q 376 388 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 389 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 377 390 ENDIF 378 391 END DO … … 437 450 ! Salinity update 438 451 ! entrapment during bottom growth 439 dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) ) & 440 & / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji) 452 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 441 453 ENDIF ! heat budget 442 454 END DO … … 476 488 zdq_i (ji ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 477 489 ENDIF 478 ! contribution to salt flux 479 zsfx_melt(ji) = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 490 ! clem: contribution to salt flux 491 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 492 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 480 493 ENDIF 481 494 END DO ! ji … … 528 541 ELSE ; zdhbf = dh_i_bott(ji) 529 542 ENDIF 543 zdvres = zdhbf - dh_i_bott(ji) 544 dh_i_bott(ji) = zdhbf 545 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 530 546 ! ! excessive energy is sent to lateral ablation 531 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) & 532 & * ( zdhbf - dh_i_bott(ji) ) * r1_rdtice 533 dh_i_bott(ji) = zdhbf 534 ! !since ice volume is only used for outputs, we keep it global for all categories 535 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 536 ! !new ice thickness 537 zhgnew (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 538 ! ! diagnostic ( bottom ice growth ) 539 ii = MOD( npb(ji) - 1, jpi ) + 1 540 ij = ( npb(ji) - 1 ) / jpi + 1 541 diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 542 diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 543 diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 547 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres * r1_rdtice 544 548 END DO 545 549 … … 552 556 ! Adapt the remaining energy if too much ice melts 553 557 !-------------------------------------------------- 554 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) ! =1 if ice 555 ! 0 if no more ice 556 zhgnew (ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 557 ! remaining heat 558 zdvres = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) ) 559 zdvsur = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first 560 zdvbot = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom 561 dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem 562 dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem 563 564 ! new ice thickness (clem) 565 zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 566 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 567 zhgnew(ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 568 569 ! !since ice volume is only used for outputs, we keep it global for all categories 570 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 571 572 ! remaining heat 558 573 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) 559 574 … … 569 584 ht_s_b(ji) = MAX( zzero , zhnfi ) 570 585 zqt_s(ji) = zqt_s(ji) * ht_s_b(ji) 586 ! we recompute dh_s_tot (clem) 587 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 571 588 572 589 ! Mass variations of ice and snow … … 579 596 ! 580 597 ! ! mass variation cumulated over category 581 rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s ! snow582 rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i ! ice598 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s ! snow 599 !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i ! ice 583 600 584 601 ! Remaining heat to the ocean … … 586 603 focea(ji) = - zfdt_final(ji) * r1_rdtice ! focea is in W.m-2 * dt 587 604 605 ! residual salt flux (clem) 606 !-------------------------- 607 ! surface 608 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvsur * rhoic * r1_rdtice 609 ! bottom 610 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting 611 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 612 ELSE ! growth 613 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 614 ENDIF 615 ! 616 ! diagnostic ( bottom ice growth ) 617 ii = MOD( npb(ji) - 1, jpi ) + 1 618 ij = ( npb(ji) - 1 ) / jpi + 1 619 diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 620 diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 621 diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 588 622 END DO 589 623 … … 591 625 592 626 !--------------------------- 593 ! Salt flux andheat fluxes627 ! heat fluxes 594 628 !--------------------------- 595 629 DO ji = kideb, kiut 596 630 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) ! =1 if ice 597 !598 ! Salt flux599 sfx_thd_1d(ji) = sfx_thd_1d(ji) + zihgnew * zsfx_melt(ji) &600 & - (1.0 - zihgnew) * zfmass_i (ji) * sm_i_b(ji) * r1_rdtice601 631 ! 602 632 ! Heat flux … … 646 676 dmgwi_1d (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 647 677 648 ! All snow is thrown in the ocean, and seawater is taken to replace the volume 649 rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic * ( 1. - rhosn / rhoic ) 650 rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 678 !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 679 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 651 680 652 681 ! Equivalent salt flux (1) Snow-ice formation component … … 658 687 ELSE ; zsm_snowice = sm_i_b(ji) 659 688 ENDIF 660 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice661 !662 689 ! entrapment during snow ice formation 663 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 664 isnowic = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch 665 IF( num_sal == 2 ) & 666 dsm_i_si_1d(ji) = ( zsm_snowice * dh_snowice(ji) & 667 & + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 ) & 668 & - sm_i_b(ji) ) * isnowic 690 ! clem: new salinity difference stored (to be used in limthd_ent.F90) 691 IF ( num_sal == 2 ) THEN 692 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) ) 693 ! salinity dif due to snow-ice formation 694 dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 695 ! salinity dif due to bottom growth 696 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0._wp ) THEN 697 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 698 ENDIF 699 ENDIF 669 700 670 701 ! Actualize new snow and ice thickness. … … 680 711 diag_sni_gr(ii,ij) = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 681 712 ! 713 ! salt flux 714 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 715 !-------------------------------- 716 ! Update mass fluxes (clem) 717 !-------------------------------- 718 rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic 719 rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn 720 682 721 END DO !ji 683 722 ! 684 723 CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 685 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, z sfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )724 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 686 725 CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) 687 726 CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) 727 ! 728 CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem 688 729 ! 689 730 END SUBROUTINE lim_thd_dh -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r3808 r4045 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 … … 34 35 35 36 !!---------------------------------------------------------------------- 36 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)37 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 37 38 !! $Id$ 38 39 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 156 157 DO ji = kideb , kiut 157 158 ! is there snow or not 158 isnow(ji)= INT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) )) )159 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) ) 159 160 ! surface temperature of fusion 160 161 !!gm ??? ztfs(ji) = rtt !!!???? 161 ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt162 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 162 163 ! layer thickness 163 zh_i(ji) = ht_i_b(ji) / nlay_i164 zh_s(ji) = ht_s_b(ji) / nlay_s164 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 165 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 165 166 END DO 166 167 … … 174 175 DO layer = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 175 176 DO ji = kideb , kiut 176 z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s177 z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s ) 177 178 END DO 178 179 END DO … … 180 181 DO layer = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 181 182 DO ji = kideb , kiut 182 z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i183 z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i ) 183 184 END DO 184 185 END DO … … 201 202 DO ji = kideb , kiut 202 203 ! switches 203 isnow(ji) = INT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) )) )204 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ) ) 204 205 ! hs > 0, isnow = 1 205 206 zhsu (ji) = hnzst ! threshold for the computation of i0 206 207 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) ) 207 208 208 i0(ji) = ( 1._wp- isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )209 i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 209 210 !fr1_i0_1d = i0 for a thin ice surface 210 211 !fr1_i0_2d = i0 for a thick ice surface … … 243 244 244 245 DO ji = kideb, kiut ! ice initialization 245 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp- isnow(ji) )246 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1 - isnow(ji) ) 246 247 END DO 247 248 … … 256 257 257 258 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)259 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 260 END DO 260 261 … … 264 265 ii = MOD( npb(ji) - 1 , jpi ) + 1 265 266 ij = ( npb(ji) - 1 ) / jpi + 1 266 fstroc(ii,ij,jl) = zradtr_i(ji,nlay_i)267 fstroc(ii,ij,jl) = iatte_1d(ji) * zradtr_i(ji,nlay_i) ! clem modif 267 268 END DO 268 269 ! +++++ … … 376 377 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & 377 378 (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 378 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)* isnow(ji) &379 + zkappa_i(ji,0)* (1.0-isnow(ji))379 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 380 + zkappa_i(ji,0)*REAL( 1 - isnow(ji) ) 380 381 END DO 381 382 ! … … 658 659 t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 659 660 * t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) & 660 * MAX(0.0,SIGN(1.0,ht_s_b(ji) -zeps))661 * MAX(0.0,SIGN(1.0,ht_s_b(ji))) 661 662 662 663 ! surface temperature 663 isnow(ji) = INT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) ) ) )664 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) ) ) ) 664 665 ztsuoldit(ji) = t_su_b(ji) 665 IF( t_su_b(ji) < ztfs(ji) ) 666 t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( isnow(ji)*t_s_b(ji,1) &667 & + (1.0-isnow(ji))*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))666 IF( t_su_b(ji) < ztfs(ji) ) & 667 t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) & 668 & + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 668 669 END DO 669 670 ! … … 721 722 #endif 722 723 ! ! surface ice conduction flux 723 isnow(ji) = INT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) )724 fc_su(ji) = - isnow(ji)* zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) &725 & - ( 1._wp- isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji))724 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) ) 725 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) & 726 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji)) 726 727 ! ! bottom ice conduction flux 727 728 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) … … 734 735 DO ji = kideb, kiut 735 736 ! Upper snow value 736 fc_s(ji,0) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) )737 fc_s(ji,0) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) ) 737 738 ! Bott. snow value 738 fc_s(ji,1) = - isnow(ji)* zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) )739 fc_s(ji,1) = - REAL( isnow(ji) ) * zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) ) 739 740 END DO 740 741 DO ji = kideb, kiut ! Upper ice layer 741 fc_i(ji,0) = - isnow(ji) * & ! interface flux if there is snow742 fc_i(ji,0) = - REAL( isnow(ji) ) * & ! interface flux if there is snow 742 743 ( zkappa_i(ji,0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) & 743 - ( 1.0- isnow(ji) ) * ( zkappa_i(ji,0) * &744 - REAL( 1 - isnow(ji) ) * ( zkappa_i(ji,0) * & 744 745 zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not 745 746 END DO -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r3625 r4045 44 44 45 45 !!---------------------------------------------------------------------- 46 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)46 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 47 47 !! $Id$ 48 48 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 75 75 76 76 INTEGER :: ji,jk ! dummy loop indices 77 INTEGER :: zji, zjj , & ! dummy indices77 INTEGER :: ii, ij , & ! dummy indices 78 78 ntop0 , & ! old layer top index 79 79 nbot1 , & ! new layer bottom index … … 145 145 146 146 DO ji = kideb, kiut 147 zh_i(ji) = old_ht_i_b(ji) / nlay_i148 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 ) 149 149 END DO 150 150 … … 166 166 DO jk = 1, nlays0 167 167 DO ji = kideb, kiut 168 snind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))) &169 + 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))))) 170 170 zdeltah(ji)= zdeltah(ji) + zh_s(ji) 171 171 END DO ! ji … … 175 175 ! 0 if not 176 176 DO ji = kideb, kiut 177 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))))) 178 178 END DO ! ji 179 179 … … 190 190 DO jk = 1, nlayi0 191 191 DO ji = kideb, kiut 192 icsuind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))) &193 + 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))))) 194 194 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 195 195 END DO ! ji … … 200 200 ! 0 if not 201 201 DO ji = kideb, kiut 202 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)) ) ) ) 203 203 ENDDO 204 204 … … 216 216 DO jk = nlayi0, 1, -1 217 217 DO ji = kideb, kiut 218 icboind(ji) = (nlayi0+1-jk) * INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))) &219 & + 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))))) 220 220 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 221 221 END DO … … 232 232 ! 0 if ablation is on the way 233 233 DO ji = kideb, kiut 234 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))))) 235 235 END DO 236 236 … … 248 248 DO ji = kideb, kiut 249 249 snicind(ji) = (nlays0+1-jk) & 250 * INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))) + snicind(ji) &251 * (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))))) 252 252 zdeltah(ji) = zdeltah(ji) + zh_s(ji) 253 253 END DO … … 258 258 ! 0 if not 259 259 DO ji = kideb, kiut 260 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))))) 261 261 ENDDO 262 262 … … 279 279 280 280 DO ji = kideb, kiut 281 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1 .- snicind(ji) ) * snicswi(ji)281 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1 - snicind(ji) ) * snicswi(ji) 282 282 ! cotes of the top of the layers 283 283 zm0(ji,0) = 0._wp … … 291 291 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 292 292 limsum = MIN( limsum , nlay_s ) 293 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * limsum294 END DO 295 END DO 296 297 DO ji = kideb, kiut 298 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - snicswi(ji) * dh_snowice(ji) + zh_s(ji) * nlays0299 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) 300 300 END DO 301 301 … … 309 309 310 310 DO ji = kideb, kiut ! layer heat content 311 qm0 (ji,1) = rhosn * ( cpic * ( rtt - ( 1.- snswi(ji) ) * tatm_ice_1d(ji) &312 & - 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) ) & 313 313 & + lfus ) * zthick0(ji,1) 314 314 zqts_in(ji) = zqts_in(ji) + qm0(ji,1) … … 320 320 limsum = MIN( limsum , nlay_s ) 321 321 qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk) 322 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, epsi20- ht_s_b(ji) ) )323 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 324 324 END DO ! jk 325 325 END DO ! ji … … 360 360 !------------------- 361 361 DO ji = kideb, kiut 362 zh_s(ji) = ht_s_b(ji) / nlay_s362 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 363 363 z_s(ji,0) = 0._wp 364 364 ENDDO … … 366 366 DO jk = 1, nlay_s 367 367 DO ji = kideb, kiut 368 z_s(ji,jk) = zh_s(ji) * jk368 z_s(ji,jk) = zh_s(ji) * REAL( jk ) 369 369 END DO 370 370 END DO … … 394 394 & - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10)) 395 395 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) & 396 & * MAX(0.0,SIGN(1.0, nbot0(ji)-layer0+epsi20))396 & * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 397 397 END DO 398 398 END DO … … 410 410 DO ji = kideb, kiut 411 411 IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN 412 zji = MOD( npb(ji) - 1, jpi ) + 1413 zjj = ( npb(ji) - 1 ) / jpi + 1412 ii = MOD( npb(ji) - 1, jpi ) + 1 413 ij = ( npb(ji) - 1 ) / jpi + 1 414 414 WRITE(numout,*) ' violation of heat conservation : ', & 415 415 ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 416 WRITE(numout,*) ' ji, jj : ', zji, zjj416 WRITE(numout,*) ' ji, jj : ', ii, ij 417 417 WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) 418 418 WRITE(numout,*) ' zqts_in : ', zqts_in (ji) * r1_rdtice … … 441 441 DO jk = 1, nlay_s 442 442 DO ji = kideb, kiut 443 zswitch = MAX ( 0.0 , SIGN ( 1.0, epsi20- ht_s_b(ji) ) )443 zswitch = MAX ( 0.0 , SIGN ( 1.0, - ht_s_b(ji) ) ) 444 444 t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 ) 445 445 END DO … … 480 480 limsum = ( (icsuswi(ji)*(icsuind(ji)+jk-1) + & 481 481 (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) 482 zm0(ji,jk)= icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) &483 + limsum* zh_i(ji)484 END DO 485 END DO 486 487 DO ji = kideb, kiut 488 zm0(ji,nbot0(ji)) = icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) + dh_i_bott(ji) &489 + zh_i(ji) * nlayi0490 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) 491 491 END DO 492 492 … … 521 521 !---------------------------- 522 522 DO ji = kideb, kiut 523 ztmelts = ( 1.0- icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice524 & + 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 525 525 & + rtt ! in Kelvin 526 526 … … 528 528 ztform = t_i_b(ji,nlay_i) 529 529 IF( num_sal == 2 ) ztform = t_bo_b(ji) 530 qm0(ji,nbot0(ji)) = ( 1.0- icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice531 & + 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 532 532 + lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) ) & 533 533 - rcp*(ztmelts-rtt) ) * zthick0(ji,nbot0(ji) ) … … 540 540 ! energy of the flooding seawater 541 541 zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * & 542 (rhoic - rhosn) / rhoic * snicswi(ji) ! generally positive542 (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive 543 543 ! Heat conservation diagnostic 544 544 qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic … … 549 549 ! = enthalpy of snow + enthalpy of frozen water 550 550 zqsnic = zqsnow(ji) + zqsnic 551 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 552 553 553 END DO ! ji … … 556 556 DO ji = kideb, kiut 557 557 ! Heat conservation 558 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06 +epsi20) ) &559 & * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + epsi20) )558 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06) ) & 559 & * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) 560 560 END DO 561 561 END DO … … 575 575 !------------------ 576 576 DO ji = kideb, kiut 577 zh_i(ji) = ht_i_b(ji) / nlay_i577 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 578 578 ENDDO 579 579 … … 606 606 q_i_b(ji,layer1) = q_i_b(ji,layer1) & 607 607 + zrl01(layer1,layer0)*qm0(ji,layer0) & 608 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06 +epsi20)) &609 * MAX(0.0,SIGN(1.0, nbot0(ji)-layer0+epsi20))608 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06)) & 609 * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 610 610 END DO 611 611 END DO … … 624 624 DO ji = kideb, kiut 625 625 IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN 626 zji = MOD( npb(ji) - 1, jpi ) + 1627 zjj = ( npb(ji) - 1 ) / jpi + 1626 ii = MOD( npb(ji) - 1, jpi ) + 1 627 ij = ( npb(ji) - 1 ) / jpi + 1 628 628 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 629 WRITE(numout,*) ' ji, jj : ', zji, zjj629 WRITE(numout,*) ' ji, jj : ', ii, ij 630 630 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 631 631 WRITE(numout,*) ' zqti_in : ', zqti_in (ji) * r1_rdtice -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r3625 r4045 46 46 47 47 !!---------------------------------------------------------------------- 48 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)48 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 49 49 !! $Id$ 50 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 78 78 !! update ht_s_b, ht_i_b and tbif_1d(:,:) 79 79 !!------------------------------------------------------------------------ 80 INTEGER 81 INTEGER 82 INTEGER :: zji, zjj, iter ! - -83 REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde! local scalars80 INTEGER :: ji,jj,jk,jl,jm ! dummy loop indices 81 INTEGER :: layer, nbpac ! local integers 82 INTEGER :: ii, ij, iter ! - - 83 REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde ! local scalars 84 84 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 85 85 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 86 REAL(wp) :: zcoef ! - -87 86 LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness 88 87 CHARACTER (len = 15) :: fieldid … … 160 159 DO ji = 1, jpi 161 160 !Energy of melting q(S,T) [J.m-3] 162 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * nlay_i161 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * REAL( nlay_i ) 163 162 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) ) ) !0 if no ice and 1 if yes 164 163 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb … … 342 341 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 343 342 DO ji = 1, nbpac 344 zji = MOD( npac(ji) - 1 , jpi ) + 1345 zjj = ( npac(ji) - 1 ) / jpi + 1346 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m( zji,zjj) )343 ii = MOD( npac(ji) - 1 , jpi ) + 1 344 ij = ( npac(ji) - 1 ) / jpi + 1 345 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij) ) 347 346 END DO 348 347 CASE ( 3 ) ! Sice = F(z) [multiyear ice] … … 389 388 END DO 390 389 391 !---------------------------------392 ! Salt flux due to new ice growth393 !---------------------------------394 ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine)395 DO ji = 1, nbpac396 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice397 rdm_ice_1d(ji) = rdm_ice_1d(ji) + rhoic * zv_newice(ji)398 END DO ! ji399 400 390 !------------------------------------ 401 391 ! Diags for energy conservation test 402 392 !------------------------------------ 403 393 DO ji = 1, nbpac 404 zji = MOD( npac(ji) - 1 , jpi ) + 1405 zjj = ( npac(ji) - 1 ) / jpi + 1394 ii = MOD( npac(ji) - 1 , jpi ) + 1 395 ij = ( npac(ji) - 1 ) / jpi + 1 406 396 ! 407 zde = ze_newice(ji) / unit_fac * area( zji,zjj) * zv_newice(ji)397 zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji) 408 398 ! 409 vt_i_init( zji,zjj) = vt_i_init(zji,zjj) + zv_newice(ji) ! volume410 et_i_init( zji,zjj) = et_i_init(zji,zjj) + zde ! Energy399 vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji) ! volume 400 et_i_init(ii,ij) = et_i_init(ii,ij) + zde ! Energy 411 401 412 402 END DO … … 419 409 !----------------- 420 410 DO ji = 1, nbpac 421 zji = MOD( npac(ji) - 1 , jpi ) + 1422 zjj = ( npac(ji) - 1 ) / jpi + 1411 ii = MOD( npac(ji) - 1 , jpi ) + 1 412 ij = ( npac(ji) - 1 ) / jpi + 1 423 413 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 424 diag_lat_gr( zji,zjj) = zv_newice(ji) * r1_rdtice414 diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem 425 415 END DO !ji 426 416 … … 441 431 ! we keep the excessive volume in memory and attribute it later to bottom accretion 442 432 DO ji = 1, nbpac 443 IF ( za_newice(ji) > ( 1._wp- zat_i_ac(ji) ) ) THEN444 zda_res(ji) = za_newice(ji) - ( 1.0- zat_i_ac(ji) )433 IF ( za_newice(ji) > ( amax - zat_i_ac(ji) ) ) THEN 434 zda_res(ji) = za_newice(ji) - ( amax - zat_i_ac(ji) ) 445 435 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 446 436 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 473 463 DO ji = 1, nbpac 474 464 jl = zcatac(ji) ! categroy in which new ice is put 475 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) ) ) ! zindb=1 if ice =0 otherwise465 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 476 466 zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness 477 467 zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness 478 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi1 1) ) ! ice totally new in jl category468 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! ice totally new in jl category 479 469 END DO 480 470 … … 482 472 DO ji = 1, nbpac 483 473 jl = zcatac(ji) 484 zqold = ze_i_ac(ji,jk,jl) 485 zalphai = MIN( zhice_old(ji,jl) * jk / nlay_i, zh_newice(ji) ) &486 & - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i, zh_newice(ji) )474 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 475 zalphai = MIN( zhice_old(ji,jl) * REAL( jk ) / REAL( nlay_i ), zh_newice(ji) ) & 476 & - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) ) 487 477 ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) & 488 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / nlay_i&478 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / REAL( nlay_i ) & 489 479 + za_newice(ji) * ze_newice(ji) * zalphai & 490 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( zv_i_ac(ji,jl) / nlay_i)480 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) ) 491 481 END DO 492 482 END DO … … 513 503 DO ji = 1, nbpac 514 504 zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) ) 515 zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 ) 505 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi06 ) ) ! clem 506 zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 ) 516 507 END DO 517 508 END DO … … 524 515 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 525 516 DO ji = 1, nbpac 526 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) ) ) ! zindb=1 if ice =0 otherwise517 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 527 518 zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 528 519 zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb & … … 536 527 DO jk = 1, nlay_i 537 528 DO ji = 1, nbpac 538 zthick0(ji,jk,jl) = zhice_old(ji,jl) / nlay_i529 zthick0(ji,jk,jl) = zhice_old(ji,jl) / REAL( nlay_i ) 539 530 zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl) 540 531 END DO … … 555 546 DO layer = 1, nlay_i + 1 556 547 DO ji = 1, nbpac 557 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) )548 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 558 549 ! Redistributing energy on the new grid 559 zweight = MAX ( MIN( zhice_old(ji,jl) * layer , zdummy(ji,jl) * jk) &560 & - MAX( zhice_old(ji,jl) * ( layer - 1 ) , zdummy(ji,jl) *( jk - 1 ) ) , 0._wp ) &561 & /( MAX( nlay_i* zthick0(ji,layer,jl),epsi10) ) * zindb550 zweight = MAX ( MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) ) & 551 & - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp ) & 552 & /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb 562 553 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 563 554 END DO ! ji … … 569 560 DO jk = 1, nlay_i 570 561 DO ji = 1, nbpac 571 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )562 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 572 563 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) & 573 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * nlay_i* zindb564 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb 574 565 END DO 575 566 END DO … … 581 572 DO jl = 1, jpl 582 573 DO ji = 1, nbpac 583 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes574 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 584 575 zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 585 576 END DO … … 589 580 ! Update salinity 590 581 !----------------- 591 IF( num_sal == 2 ) THEN ! Sice = F(z,t)582 !clem IF( num_sal == 2 ) THEN 592 583 DO jl = 1, jpl 593 584 DO ji = 1, nbpac 594 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )! 0 if no ice and 1 if yes585 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 595 586 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 596 zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb587 zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif 597 588 END DO 598 589 END DO 599 ENDIF 590 !clem ENDIF 591 592 !-------------------------------- 593 ! Update mass/salt fluxes (clem) 594 !-------------------------------- 595 DO jl = 1, jpl 596 DO ji = 1, nbpac 597 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 598 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 599 rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb 600 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb 601 END DO 602 END DO 600 603 601 604 !------------------------------------------------------------------------------! … … 606 609 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 607 610 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 608 IF ( num_sal == 2 ) &611 !clem IF ( num_sal == 2 ) & 609 612 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 610 613 DO jk = 1, nlay_i … … 622 625 DO jl = 1, jpl 623 626 DO jk = 1, nlay_i ! heat content in 10^9 Joules 624 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i/ unit_fac627 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac 625 628 END DO 626 629 END DO -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
r3625 r4045 33 33 34 34 !!---------------------------------------------------------------------- 35 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)35 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 36 36 !! $Id$ 37 37 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 82 82 DO ji = kideb, kiut 83 83 zhiold(ji) = ht_i_b(ji) - dh_i_bott(ji) - dh_snowice(ji) - dh_i_surf(ji) 84 zsiold(ji) = sm_i_b(ji) 84 85 END DO 85 86 … … 90 91 DO jk = 1, nlay_i 91 92 DO ji = kideb, kiut 92 ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i93 ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / REAL (nlay_i ) 93 94 END DO 94 95 END DO … … 117 118 ! only drainage terms ( gravity drainage and flushing ) 118 119 ! snow ice / bottom sources are added in lim_thd_ent to conserve energy 119 zsiold(ji) = sm_i_b(ji)120 120 sm_i_b(ji) = sm_i_b(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji) 121 121 … … 123 123 i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) ) 124 124 sm_i_b(ji) = i_ice_switch * sm_i_b(ji) + s_i_min * ( 1._wp - i_ice_switch ) 125 END DO ! ji 126 127 CALL lim_var_salprof1d( kideb, kiut ) ! Salinity profile 128 129 130 !---------------------------- 131 ! Heat flux - brine drainage 132 !---------------------------- 133 134 DO ji = kideb, kiut 135 !!gm useless 136 ! iflush : 1 if summer 137 iflush = MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt ) ) 138 ! igravdr : 1 if t_su lt t_bo 139 igravdr = MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) ) ) 140 ! iaccrbo : 1 if bottom accretion 141 iaccrbo = MAX( 0._wp , SIGN( 1._wp , dh_i_bott(ji) ) ) 142 !!gm end useless 143 ! 125 126 !---------------------------- 127 ! Heat flux - brine drainage 128 !---------------------------- 144 129 fhbri_1d(ji) = 0._wp 145 END DO ! ji 146 147 !---------------------------- 148 ! Salt flux - brine drainage 149 !---------------------------- 150 DO ji = kideb, kiut 151 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) ) 152 sfx_bri_1d(ji) = sfx_bri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) & 153 & * ( MAX( dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji) , sm_i_b(ji) - zsiold(ji) ) ) * r1_rdtice 154 END DO 130 131 !---------------------------- 132 ! Salt flux - brine drainage 133 !---------------------------- 134 sfx_bri_1d(ji) = sfx_bri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) * r1_rdtice 135 136 END DO 137 138 ! Salinity profile 139 CALL lim_var_salprof1d( kideb, kiut ) 140 155 141 156 142 ! Only necessary for conservation check since salinity is modified … … 178 164 IF( num_sal == 3 ) CALL lim_var_salprof1d( kideb, kiut ) 179 165 180 181 !------------------------------------------------------------------------------|182 ! 5) Computation of salt flux due to Bottom growth183 !------------------------------------------------------------------------------|184 ! note: s_i_new = bulk_sal in constant salinity case185 DO ji = kideb, kiut186 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0._wp ) * r1_rdtice187 END DO188 166 ! 189 167 CALL wrk_dealloc( jpij, ze_init, zhiold, zsiold ) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r3625 r4045 28 28 USE prtctl ! Print control 29 29 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 30 USE limvar ! clem for ice thickness correction 30 31 31 32 IMPLICIT NONE … … 36 37 REAL(wp) :: epsi06 = 1.e-06_wp ! constant values 37 38 REAL(wp) :: epsi03 = 1.e-03_wp 38 REAL(wp) :: zeps10 = 1.e-10_wp39 REAL(wp) :: epsi10 = 1.e-10_wp 39 40 REAL(wp) :: epsi16 = 1.e-16_wp 41 REAL(wp) :: epsi20 = 1.e-20_wp 40 42 REAL(wp) :: rzero = 0._wp 41 43 REAL(wp) :: rone = 1._wp … … 46 48 # include "vectopt_loop_substitute.h90" 47 49 !!---------------------------------------------------------------------- 48 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)50 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 49 51 !! $Id$ 50 52 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 69 71 INTEGER :: initad ! number of sub-timestep for the advection 70 72 INTEGER :: ierr ! error status 71 REAL(wp) :: zindb , zindsn , zindic ! local scalar73 REAL(wp) :: zindb , zindsn , zindic, zindh, zinda ! local scalar 72 74 REAL(wp) :: zusvosn, zusvoic, zbigval ! - - 73 75 REAL(wp) :: zcfl , zusnit , zrtt ! - - … … 77 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 78 80 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e 81 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) 82 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (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 and thickness 87 REAL(wp) :: zdv, zda, zvi, zvs, zsmv 79 88 !!--------------------------------------------------------------------- 80 89 … … 82 91 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 83 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( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 103 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 104 ENDIF 105 !- check conservation (C Rousset) 106 ! ------------------------------- 84 107 85 108 IF( numit == nstart .AND. lwp ) THEN … … 96 119 IF( ln_limdyn ) THEN ! Advection of sea ice properties ! 97 120 ! !-------------------------------------! 98 ! 99 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 !zhimax(ji,jj,jl) = ( ht_i(ji ,jj ,jl) * tmask(ji, jj ,1) + ht_i(ji-1,jj-1,jl) * tmask(ji-1,jj-1,1) + ht_i(ji+1,jj+1,jl) * tmask(ji+1,jj+1,1) & 137 ! & + ht_i(ji-1,jj ,jl) * tmask(ji-1,jj ,1) + ht_i(ji ,jj-1,jl) * tmask(ji ,jj-1,1) & 138 ! & + ht_i(ji+1,jj ,jl) * tmask(ji+1,jj ,1) + ht_i(ji ,jj+1,jl) * tmask(ji ,jj+1,1) & 139 ! & + ht_i(ji-1,jj+1,jl) * tmask(ji-1,jj+1,1) + ht_i(ji+1,jj-1,jl) * tmask(ji+1,jj-1,1) ) 140 END DO 141 END DO 142 CALL lbc_lnk(zhimax(:,:,jl),'T',1.) 143 END DO 144 100 145 !------------------------- 101 146 ! transported fields … … 126 171 ! ENDIF 127 172 !!gm end 128 initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )173 initad = 1 + NINT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) ) 129 174 zusnit = 1.0 / REAL( initad ) 130 175 IF( zcfl > 0.5 .AND. lwp ) & … … 282 327 END DO 283 328 284 !-----------------------------------------285 ! Remultiply everything by ice area286 !-----------------------------------------287 zs0ow(:,:) = MAX( rzero, zs0ow(:,:) * area(:,:) )288 DO jl = 1, jpl289 zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) ) !!bug: est-ce utile290 zs0sn (:,:,jl) = MAX( rzero, zs0sn (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres291 zs0sm (:,:,jl) = MAX( rzero, zs0sm (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres292 zs0oi (:,:,jl) = MAX( rzero, zs0oi (:,:,jl) * area(:,:) )293 zs0a (:,:,jl) = MAX( rzero, zs0a (:,:,jl) * area(:,:) ) !! suppress both change le resultat294 zs0c0 (:,:,jl) = MAX( rzero, zs0c0 (:,:,jl) * area(:,:) )295 DO jk = 1, nlay_i296 zs0e(:,:,jk,jl) = MAX( rzero, zs0e (:,:,jk,jl) * area(:,:) )297 END DO ! jk298 END DO ! jl299 300 329 !------------------------------------------------------------------------------! 301 330 ! 5) Update and limit ice properties after transport … … 305 334 ! 5.1) Recover mean values over the grid squares. 306 335 !-------------------------------------------------- 307 308 DO jl = 1, jpl309 DO jk = 1, nlay_i310 DO jj = 1, jpj311 DO ji = 1, jpi312 zs0e(ji,jj,jk,jl) = MAX( rzero, zs0e(ji,jj,jk,jl) / area(ji,jj) )313 END DO314 END DO315 END DO316 END DO317 318 DO jj = 1, jpj319 DO ji = 1, jpi320 zs0ow(ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) )321 END DO322 END DO323 324 336 zs0at(:,:) = 0._wp 325 337 DO jl = 1, jpl 326 338 DO jj = 1, jpj 327 339 DO ji = 1, jpi 328 zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl) /area(ji,jj))329 zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl) /area(ji,jj))330 zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl) /area(ji,jj))331 zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl) /area(ji,jj))332 zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl) /area(ji,jj))333 zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl) /area(ji,jj))340 zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl) ) 341 zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl) ) 342 zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl) ) 343 zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl) ) 344 zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl) ) 345 zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl) ) 334 346 zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) 335 347 END DO … … 342 354 DO jj = 1, jpj 343 355 DO ji = 1, jpi 344 zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - zeps10) )356 zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - epsi10) ) 345 357 zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp ) 346 358 ato_i(ji,jj) = zs0ow(ji,jj) … … 351 363 DO jj = 1, jpj 352 364 DO ji = 1, jpi 353 zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - zeps10) ) 365 zvi = zs0ice(ji,jj,jl) 366 zvs = zs0sn(ji,jj,jl) 354 367 ! 355 zs0a(ji,jj,jl) = zindb * MIN( zs0a(ji,jj,jl), 0.99 ) 368 zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - epsi10) ) 369 ! 356 370 v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl) 357 371 v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl) 358 372 ! 359 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )360 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )373 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 374 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 361 375 zindb = MAX( zindsn, zindic ) 376 ! 362 377 zs0a(ji,jj,jl) = zindb * zs0a(ji,jj,jl) !ice concentration 363 378 a_i (ji,jj,jl) = zs0a(ji,jj,jl) 364 379 v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl) 365 380 v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl) 381 ! 382 ! Update mass fluxes (clem) 383 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 384 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 385 END DO 386 END DO 387 END DO 388 389 !--- Thickness correction in case too high (clem) -------------------------------------------------------- 390 CALL lim_var_glo2eqv 391 DO jl = 1, jpl 392 DO jj = 1, jpj 393 DO ji = 1, jpi 394 395 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 396 zvi = v_i(ji,jj,jl) 397 zvs = v_s(ji,jj,jl) 398 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 399 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 400 401 zindh = 1._wp 402 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 403 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 404 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 405 zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 406 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 ) 407 ELSE 408 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 409 zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 410 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 ) 411 ENDIF 412 413 ! small correction due to *zindh for a_i 414 v_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) 415 v_s(ji,jj,jl) = zindh * v_s(ji,jj,jl) 416 417 ! Update mass fluxes 418 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 419 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 420 421 ENDIF 422 423 diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 424 366 425 END DO 367 426 END DO 368 427 END DO 369 428 429 ! --- 370 430 DO jj = 1, jpj 371 431 DO ji = 1, jpi 372 zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) 432 zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) ! clem@useless?? 373 433 END DO 374 434 END DO … … 378 438 !---------------------- 379 439 380 zbigval = 1. d+13440 zbigval = 1.e+13 381 441 382 442 DO jl = 1, jpl 383 443 DO jj = 1, jpj 384 444 DO ji = 1, jpi 445 zsmv = zs0sm(ji,jj,jl) 385 446 386 447 ! Switches and dummy variables … … 388 449 zusvoic = 1.0/MAX( v_i(ji,jj,jl) , epsi16 ) 389 450 zrtt = 173.15 * rone 390 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )391 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )451 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 452 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 392 453 zindb = MAX( zindsn, zindic ) 393 454 … … 395 456 zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj) , & 396 457 & zusvoic * zs0sm(ji,jj,jl) ) , s_i_min ) * v_i(ji,jj,jl) 397 IF( num_sal == 2 ) smv_i(ji,jj,jl) = zindic * zsal + (1.0-zindic) * 0._wp458 IF( num_sal == 2 ) smv_i(ji,jj,jl) = zindic * zsal 398 459 399 460 zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi16 ) ), 0._wp ) * a_i(ji,jj,jl) … … 402 463 ! Snow heat content 403 464 ze = MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval ) 404 e_s(ji,jj,1,jl) = zindsn * ze + (1.0 - zindsn) * 0.0 405 465 e_s(ji,jj,1,jl) = zindsn * ze 466 467 ! Update salt fluxes (clem) 468 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 406 469 END DO !ji 407 470 END DO !jj … … 413 476 DO ji = 1, jpi 414 477 ! Ice heat content 415 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )478 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 416 479 ze = MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval ) 417 e_i(ji,jj,jk,jl) = zindic * ze + ( 1.0 - zindic ) * 0.0480 e_i(ji,jj,jk,jl) = zindic * ze 418 481 END DO !ji 419 482 END DO ! jj 420 483 END DO ! jk 421 484 END DO ! jl 485 486 487 ! --- agglomerate variables (clem) ----------------- 488 vt_i (:,:) = 0._wp 489 vt_s (:,:) = 0._wp 490 at_i (:,:) = 0._wp 491 ! 492 DO jl = 1, jpl 493 DO jj = 1, jpj 494 DO ji = 1, jpi 495 ! 496 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 497 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 498 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 499 ! 500 zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi16 ) ) 501 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda ! ice thickness 502 END DO 503 END DO 504 END DO 505 ! ------------------------------------------------- 506 507 422 508 423 509 ENDIF … … 454 540 END DO 455 541 ENDIF 542 ! ------------------------------- 543 !- check conservation (C Rousset) 544 IF (ln_limdiahsb) THEN 545 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 546 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 547 548 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 549 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 550 551 zchk_vmin = glob_min(v_i) 552 zchk_amax = glob_max(SUM(a_i,dim=3)) 553 zchk_amin = glob_min(a_i) 554 555 IF(lwp) THEN 556 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limtrp) = ',(zchk_v_i * rday) 557 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * rday) 558 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limtrp) = ',(zchk_vmin * 1.e-3) 559 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limtrp) = ',zchk_amin 560 ENDIF 561 ENDIF 562 !- check conservation (C Rousset) 563 ! ------------------------------- 456 564 ! 457 565 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow ) 458 566 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 459 567 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e ) 568 569 CALL wrk_dealloc( jpi,jpj,jpl,zaiold, zhimax ) ! clem 460 570 ! 461 571 END SUBROUTINE lim_trp -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r3625 r4045 62 62 PUBLIC lim_var_eqv2glo ! 63 63 PUBLIC lim_var_salprof ! 64 PUBLIC lim_var_icetm ! 64 65 PUBLIC lim_var_bv ! 65 66 PUBLIC lim_var_salprof1d ! 66 67 67 REAL(wp) :: eps 20 = 1.e-20_wp ! module constants68 REAL(wp) :: eps 16 = 1.e-16_wp ! - -69 REAL(wp) :: eps 13 = 1.e-13_wp ! - -70 REAL(wp) :: eps 10 = 1.e-10_wp ! - -71 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 ! - - 72 73 REAL(wp) :: zzero = 0.e0 ! - - 73 74 REAL(wp) :: zone = 1.e0 ! - - 74 75 75 76 !!---------------------------------------------------------------------- 76 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)77 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 77 78 !! $Id$ 78 79 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 97 98 ! 98 99 INTEGER :: ji, jj, jk, jl ! dummy loop indices 99 REAL(wp) :: zinda 100 REAL(wp) :: zinda, zindb 100 101 !!------------------------------------------------------------------ 101 102 … … 116 117 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 117 118 ! 118 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10) )119 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 120 121 END DO 121 122 END DO … … 137 138 DO jj = 1, jpj 138 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 ) ) 139 142 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 140 zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - 0.10 ) ) 141 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , eps13 ) * zinda ! ice salinity 142 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) ) 143 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , eps13 ) * zinda ! ice age 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 144 145 END DO 145 146 END DO … … 175 176 DO jj = 1, jpj 176 177 DO ji = 1, jpi 177 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) ) ) !0 if no ice and 1 if yes 178 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps10 ) * zindb 179 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps10 ) * zindb 180 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps10 ) * zindb 178 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 182 a_i(ji,jj,jl) = a_i (ji,jj,jl) * zindb ! clem correction 181 183 END DO 182 184 END DO … … 187 189 DO jj = 1, jpj 188 190 DO ji = 1, jpi 189 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) ) ) !0 if no ice and 1 if yes190 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , eps 10 ) * zindb191 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 192 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * zindb 191 193 END DO 192 194 END DO … … 208 210 DO ji = 1, jpi 209 211 ! ! Energy of melting q(S,T) [J.m-3] 210 zq_i = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , eps 06 ) * REAL(nlay_i,wp)212 zq_i = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi06 ) * REAL(nlay_i,wp) 211 213 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) ) ! zindb = 0 if no ice and 1 if yes 212 214 zq_i = zq_i * unit_fac * zindb !convert units … … 234 236 DO ji = 1, jpi 235 237 !Energy of melting q(S,T) [J.m-3] 236 zq_s = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , eps 06 ) ) * REAL(nlay_s,wp)238 zq_s = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL(nlay_s,wp) 237 239 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) ) ) ! zindb = 0 if no ice and 1 if yes 238 240 zq_s = zq_s * unit_fac * zindb ! convert units … … 253 255 DO jj = 1, jpj 254 256 DO ji = 1, jpi 255 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , -a_i(ji,jj,jl) ) ) ) & 256 & * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) ) ) ) 257 tm_i(ji,jj) = tm_i(ji,jj) + t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 258 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , eps10 ) ) 257 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 258 tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 259 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 259 260 END DO 260 261 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 … … 429 458 ! 430 459 INTEGER :: ji, jk ! dummy loop indices 431 INTEGER :: zji, zjj ! local integers460 INTEGER :: ii, ij ! local integers 432 461 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal ! local scalars 433 462 REAL(wp) :: zalpha, zind0, zind01, zindbal, zs_zero ! - - … … 463 492 !CDIR NOVERRCHK 464 493 DO ji = kideb, kiut 465 zji = MOD( npb(ji) - 1 , jpi ) + 1466 zjj = ( npb(ji) - 1 ) / jpi + 1494 ii = MOD( npb(ji) - 1 , jpi ) + 1 495 ij = ( npb(ji) - 1 ) / jpi + 1 467 496 ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise 468 497 zind0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i_b(ji) ) ) … … 470 499 zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) ) 471 500 ! if 2.sm_i GE sss_m then zindbal = 1 472 zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m( zji,zjj) ) )501 zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(ii,ij) ) ) 473 502 ! 474 503 zalpha = ( zind0 + zind01 * ( sm_i_b(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zindbal ) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r3853 r4045 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 27 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 28 USE iom 29 USE timing ! Timing 30 USE lib_fortran ! Fortran utilities 28 31 29 32 IMPLICIT NONE … … 31 34 32 35 PUBLIC lim_wri ! routine called by lim_step.F90 33 34 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 35 39 36 40 INTEGER :: noumef ! number of fields … … 48 52 INTEGER , DIMENSION(jpnoumax) :: nc , nca ! switch for saving field ( = 1 ) or not ( = 0 ) 49 53 50 REAL(wp) :: epsi 16 = 1e-16_wp54 REAL(wp) :: epsi06 = 1e-6_wp 51 55 REAL(wp) :: zzero = 0._wp 52 56 REAL(wp) :: zone = 1._wp 53 57 !!---------------------------------------------------------------------- 54 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)58 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 55 59 !! $Id$ 56 60 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 77 81 INTEGER :: ierr 78 82 REAL(wp),DIMENSION(1) :: zdept 79 REAL(wp) :: zsto, zjulian, zout, zindh, zinda, zindb 83 REAL(wp) :: zsto, zjulian, zout, zindh, zinda, zindb, zindc 80 84 REAL(wp), POINTER, DIMENSION(:,:,:) :: zcmo, zcmoa 81 85 REAL(wp), POINTER, DIMENSION(:,: ) :: zfield 82 86 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmaskitd, zoi, zei 83 87 84 CHARACTER(len = 40) :: clhstnam, clop, clhstnama88 CHARACTER(len = 60) :: clhstnam, clop, clhstnama 85 89 86 90 INTEGER , SAVE :: nice, nhorid, ndim, niter, ndepid … … 90 94 !!------------------------------------------------------------------- 91 95 96 IF( nn_timing == 1 ) CALL timing_start('limwri') 97 92 98 CALL wrk_alloc( jpi, jpj, zfield ) 93 99 CALL wrk_alloc( jpi, jpj, jpnoumax, zcmo, zcmoa ) … … 116 122 ! Normal file 117 123 !------------- 118 119 zsto = rdt_ice120 IF( ln_mskland ) THEN ; clop = "ave(only(x))" ! put 1.e+20 on land (very expensive!!)121 ELSE ; clop = "ave(x)" ! no use of the mask value (require less cpu time)122 ENDIF123 zout = nwrite * rdt_ice / nn_fsbc124 124 niter = ( nit000 - 1 ) / nn_fsbc 125 zdept(1) = 0.126 127 125 CALL ymds2ju ( nyear, nmonth, nday, rdt, zjulian ) 128 126 zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment 129 CALL dia_nam ( clhstnam, nwrite, 'icemod' ) 130 CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit, 1, jpi, 1, jpj, niter, zjulian, rdt_ice, & 131 & nhorid, nice, domain_id=nidom, snc4chunks=snc4set ) 132 CALL histvert( nice, "deptht", "Vertical T levels", "m", 1, zdept, ndepid, "down") 133 CALL wheneq ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim) 134 135 DO jf = 1 , noumef 136 IF(lwp) WRITE(numout,*) 'jf', jf 137 IF ( nc(jf) == 1 ) THEN 138 CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj & 139 , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout ) 140 IF(lwp) WRITE(numout,*) 'nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout' 141 IF(lwp) WRITE(numout,*) nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout 142 ENDIF 143 END DO 144 145 CALL histend(nice, snc4set) 146 127 !clem 128 ! zsto = rdt_ice 129 ! IF( ln_mskland ) THEN ; clop = "ave(only(x))" ! put 1.e+20 on land (very expensive!!) 130 ! ELSE ; clop = "ave(x)" ! no use of the mask value (require less cpu time) 131 ! ENDIF 132 ! zout = nwrite * rdt_ice / nn_fsbc 133 ! zdept(1) = 0. 134 ! 135 ! CALL dia_nam ( clhstnam, nwrite, 'icemod_old' ) 136 ! CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit, 1, jpi, 1, jpj, niter, zjulian, rdt_ice, & 137 ! & nhorid, nice, domain_id=nidom, snc4chunks=snc4set ) 138 ! CALL histvert( nice, "deptht", "Vertical T levels", "m", 1, zdept, ndepid, "down") 139 ! CALL wheneq ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim) 140 ! 141 ! DO jf = 1 , noumef 142 ! IF(lwp) WRITE(numout,*) 'jf', jf 143 ! IF ( nc(jf) == 1 ) THEN 144 ! CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj & 145 ! , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout ) 146 ! IF(lwp) WRITE(numout,*) 'nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout' 147 ! IF(lwp) WRITE(numout,*) nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout 148 ! ENDIF 149 ! END DO 150 ! 151 ! CALL histend(nice, snc4set) 152 !clem 153 ! 147 154 !----------------- 148 155 ! ITD file output … … 159 166 nhorida, & ! ? linked with horizontal ... 160 167 nicea , domain_id=nidom, snc4chunks=snc4set) ! file 161 CALL histvert( nicea, "icethi", "L levels", & 162 "m", ipl , hi_mean , nz ) 168 CALL histvert( nicea, "icethi", "L levels","m", ipl , hi_mean , nz ) 163 169 DO jl = 1, jpl 164 170 zmaskitd(:,:,jl) = tmask(:,:,1) … … 198 204 zcmoa( 1:jpi, 1:jpj, 1:jpnoumax ) = 0._wp 199 205 206 ! Ice surface temperature and some fluxes 200 207 DO jl = 1, jpl 201 208 DO jj = 1, jpj 202 209 DO ji = 1, jpi 203 zindh = MAX( zzero , SIGN( zone , vt_i(ji,jj) * at_i(ji,jj) - 0.10 ) ) 204 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) ) 210 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi06 ) ) 205 211 zcmo(ji,jj,17) = zcmo(ji,jj,17) + a_i(ji,jj,jl)*qsr_ice (ji,jj,jl) 206 212 zcmo(ji,jj,18) = zcmo(ji,jj,18) + a_i(ji,jj,jl)*qns_ice(ji,jj,jl) 207 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 213 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) 214 zcmo(ji,jj,21) = zcmo(ji,jj,21) + zinda*oa_i(ji,jj,jl)/MAX(at_i(ji,jj),epsi06) 208 215 END DO 209 216 END DO 210 217 END DO 211 218 219 ! Mean sea ice temperature 220 CALL lim_var_icetm 221 222 ! Brine volume 212 223 CALL lim_var_bv 213 224 214 225 DO jj = 2 , jpjm1 215 226 DO ji = 2 , jpim1 216 zindh = MAX( zzero , SIGN( zone , vt_i(ji,jj) * at_i(ji,jj) - 0.10 ) ) 217 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) ) 218 zindb = zindh * zinda 227 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi06 ) ) 228 zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) ) ) 219 229 220 230 zcmo(ji,jj,1) = at_i(ji,jj) 221 zcmo(ji,jj,2) = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi 16 ) * zinda222 zcmo(ji,jj,3) = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi 16 ) * zinda223 zcmo(ji,jj,4) = diag_bot_gr(ji,jj) * 86400.0 * zinda! Bottom thermodynamic ice production224 zcmo(ji,jj,5) = diag_dyn_gr(ji,jj) * 86400.0 * zinda! Dynamic ice production (rid/raft)225 zcmo(ji,jj,22) = diag_lat_gr(ji,jj) * 86400.0 * zinda! Lateral thermodynamic ice production226 zcmo(ji,jj,23) = diag_sni_gr(ji,jj) * 86400.0 * zinda! Snow ice production ice production227 zcmo(ji,jj,24) = tm_i(ji,jj) - rtt228 229 zcmo(ji,jj,6) = fbif 230 zcmo(ji,jj,7) = zindb *( u_ice(ji,jj) * tmu(ji,jj) + u_ice(ji-1,jj) * tmu(ji-1,jj) ) * 0.5_wp231 zcmo(ji,jj,8) = zindb *( v_ice(ji,jj) * tmv(ji,jj) + v_ice(ji,jj-1) * tmv(ji,jj-1) ) * 0.5_wp231 zcmo(ji,jj,2) = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zinda 232 zcmo(ji,jj,3) = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zinda 233 zcmo(ji,jj,4) = diag_bot_gr(ji,jj) * rday ! Bottom thermodynamic ice production 234 zcmo(ji,jj,5) = diag_dyn_gr(ji,jj) * rday ! Dynamic ice production (rid/raft) 235 zcmo(ji,jj,22) = diag_lat_gr(ji,jj) * rday ! Lateral thermodynamic ice production 236 zcmo(ji,jj,23) = diag_sni_gr(ji,jj) * rday ! Snow ice production ice production 237 zcmo(ji,jj,24) = (tm_i(ji,jj) - rtt) * zinda 238 239 zcmo(ji,jj,6) = fbif(ji,jj)*at_i(ji,jj) 240 zcmo(ji,jj,7) = ( u_ice(ji,jj) * tmu(ji,jj) + u_ice(ji-1,jj) * tmu(ji-1,jj) ) * 0.5_wp 241 zcmo(ji,jj,8) = ( v_ice(ji,jj) * tmv(ji,jj) + v_ice(ji,jj-1) * tmv(ji,jj-1) ) * 0.5_wp 232 242 zcmo(ji,jj,9) = sst_m(ji,jj) 233 243 zcmo(ji,jj,10) = sss_m(ji,jj) … … 243 253 zcmo(ji,jj,19) = sprecip(ji,jj) 244 254 zcmo(ji,jj,20) = smt_i(ji,jj) 245 zcmo(ji,jj,21) = ot_i(ji,jj)246 255 zcmo(ji,jj,25) = et_i(ji,jj) 247 256 zcmo(ji,jj,26) = et_s(ji,jj) … … 250 259 251 260 zcmo(ji,jj,30) = bv_i(ji,jj) 252 zcmo(ji,jj,31) = hicol(ji,jj) 261 zcmo(ji,jj,31) = hicol(ji,jj) * zindb 253 262 zcmo(ji,jj,32) = strength(ji,jj) 254 263 zcmo(ji,jj,33) = SQRT( zcmo(ji,jj,7)*zcmo(ji,jj,7) + zcmo(ji,jj,8)*zcmo(ji,jj,8) ) 255 zcmo(ji,jj,34) = diag_sur_me(ji,jj) * 86400.0 * zinda! Surface melt256 zcmo(ji,jj,35) = diag_bot_me(ji,jj) * 86400.0 * zinda! Bottom melt264 zcmo(ji,jj,34) = diag_sur_me(ji,jj) * rday ! Surface melt 265 zcmo(ji,jj,35) = diag_bot_me(ji,jj) * rday ! Bottom melt 257 266 zcmo(ji,jj,36) = divu_i(ji,jj) 258 267 zcmo(ji,jj,37) = shear_i(ji,jj) 259 END DO 268 zcmo(ji,jj,38) = diag_res_pr(ji,jj) * rday ! Bottom melt 269 zcmo(ji,jj,39) = vt_i(ji,jj) ! ice volume 270 zcmo(ji,jj,40) = vt_s(ji,jj) ! snow volume 271 272 zcmo(ji,jj,41) = sfx_mec(ji,jj) 273 zcmo(ji,jj,42) = sfx_res(ji,jj) 274 275 zcmo(ji,jj,43) = diag_trp_vi(ji,jj) * rday ! transport of ice volume 276 277 END DO 260 278 END DO 261 279 … … 264 282 ! 265 283 niter = niter + 1 266 DO jf = 1 , noumef 267 ! 268 zfield(:,:) = zcmo(:,:,jf) * cmulti(jf) + cadd(jf) 269 ! 270 IF( jf == 7 .OR. jf == 8 .OR. jf == 15 .OR. jf == 16 ) THEN ; CALL lbc_lnk( zfield, 'T', -1. ) 271 ELSE ; CALL lbc_lnk( zfield, 'T', 1. ) 272 ENDIF 273 ! 274 IF( ln_nicep ) THEN 275 WRITE(numout,*) 276 WRITE(numout,*) 'nc(jf), nice, nam(jf), niter, ndim' 277 WRITE(numout,*) nc(jf), nice, nam(jf), niter, ndim 278 ENDIF 279 IF( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 ) 280 ! 281 END DO 282 283 IF( ( nn_fsbc * niter ) >= nitend .OR. kindic < 0 ) THEN 284 IF( lwp) WRITE(numout,*) ' Closing the icemod file ' 285 CALL histclo( nice ) 286 ENDIF 284 !clem 285 ! DO jf = 1 , noumef 286 ! ! 287 ! zfield(:,:) = zcmo(:,:,jf) * cmulti(jf) + cadd(jf) 288 ! ! 289 ! IF( jf == 7 .OR. jf == 8 .OR. jf == 15 .OR. jf == 16 ) THEN ; CALL lbc_lnk( zfield, 'T', -1. ) 290 ! ELSE ; CALL lbc_lnk( zfield, 'T', 1. ) 291 ! ENDIF 292 ! ! 293 ! IF( ln_nicep ) THEN 294 ! WRITE(numout,*) 295 ! WRITE(numout,*) 'nc(jf), nice, nam(jf), niter, ndim' 296 ! WRITE(numout,*) nc(jf), nice, nam(jf), niter, ndim 297 ! ENDIF 298 ! IF( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 ) 299 ! ! 300 ! END DO 301 ! 302 ! IF( ( nn_fsbc * niter ) >= nitend .OR. kindic < 0 ) THEN 303 ! IF( lwp) WRITE(numout,*) ' Closing the icemod file ' 304 ! CALL histclo( nice ) 305 ! ENDIF 306 !clem 307 ! 308 CALL iom_put ('iceconc', zcmo(:,:,1) ) ! field1: ice concentration 309 CALL iom_put ('icethic_cea', zcmo(:,:,2) ) ! field2: ice thickness (i.e. icethi(:,:)) 310 CALL iom_put ('snowthic_cea', zcmo(:,:,3)) ! field3: snow thickness 311 CALL iom_put ('icebopr', zcmo(:,:,4) ) ! field4: daily bottom thermo ice production 312 CALL iom_put ('icedypr', zcmo(:,:,5) ) ! field5: daily dynamic ice production 313 CALL iom_put ('ioceflxb', zcmo(:,:,6) ) ! field6: Oceanic flux at the ice base 314 CALL iom_put ('uice_ipa', zcmo(:,:,7) ) ! field7: ice velocity u component 315 CALL iom_put ('vice_ipa', zcmo(:,:,8) ) ! field8: ice velocity v component 316 CALL iom_put ('isst', zcmo(:,:,9) ) ! field 9: sea surface temperature 317 CALL iom_put ('isss', zcmo(:,:,10) ) ! field 10: sea surface salinity 318 CALL iom_put ('qt_oce', zcmo(:,:,11) ) ! field 11: total flux at ocean surface 319 CALL iom_put ('qsr_oce', zcmo(:,:,12) ) ! field 12: solar flux at ocean surface 320 CALL iom_put ('qns_oce', zcmo(:,:,13) ) ! field 13: non-solar flux at ocean surface 321 !CALL iom_put ('hfbri', fhbri ) ! field 14: heat flux due to brine release 322 CALL iom_put( 'utau_ice', zcmo(:,:,15) ) ! Wind stress over ice along i-axis at I-point 323 CALL iom_put( 'vtau_ice', zcmo(:,:,16) ) ! Wind stress over ice along j-axis at I-point 324 CALL iom_put ('qsr_io', zcmo(:,:,17) ) ! field 17: solar flux at ice/ocean surface 325 CALL iom_put ('qns_io', zcmo(:,:,18) ) ! field 18: non-solar flux at ice/ocean surface 326 !CALL iom_put ('snowpre', zcmo(:,:,19) * rday ! field 19 :snow precip 327 CALL iom_put ('micesalt', zcmo(:,:,20) ) ! field 20 :mean ice salinity 328 CALL iom_put ('miceage', zcmo(:,:,21) / 365) ! field 21: mean ice age 329 CALL iom_put ('icelapr',zcmo(:,:,22) ) ! field 22: daily lateral thermo ice prod. 330 CALL iom_put ('icesipr',zcmo(:,:,23) ) ! field 23: daily snowice ice prod. 331 CALL iom_put ('micet', zcmo(:,:,24) ) ! field 24: mean ice temperature 332 CALL iom_put ('icehc', zcmo(:,:,25) ) ! field 25: ice total heat content 333 CALL iom_put ('isnowhc', zcmo(:,:,26) ) ! field 26: snow total heat content 334 CALL iom_put ('icest', zcmo(:,:,27) ) ! field 27: ice surface temperature 335 CALL iom_put ('sfxbri', zcmo(:,:,28) * rday ) ! field 28: brine salt flux 336 CALL iom_put ('sfxthd', zcmo(:,:,29) * rday ) ! field 29: equivalent FW salt flux 337 CALL iom_put ('ibrinv', zcmo(:,:,30) *100 ) ! field 30: brine volume 338 CALL iom_put ('icecolf', zcmo(:,:,31) ) ! field 31: frazil ice collection thickness 339 CALL iom_put ('icestr', zcmo(:,:,32) * 0.001 ) ! field 32: ice strength 340 CALL iom_put ('icevel', zcmo(:,:,33) ) ! field 33: ice velocity 341 CALL iom_put ('isume', zcmo(:,:,34) ) ! field 34: surface melt 342 CALL iom_put ('ibome', zcmo(:,:,35) ) ! field 35: bottom melt 343 CALL iom_put ('idive', zcmo(:,:,36) * 1.0e8) ! field 36: divergence 344 CALL iom_put ('ishear', zcmo(:,:,37) * 1.0e8 ) ! field 37: shear 345 CALL iom_put ('icerepr', zcmo(:,:,38) ) ! field 38: daily prod./melting due to limupdate 346 CALL iom_put ('icevolu', zcmo(:,:,39) ) ! field 39: ice volume 347 CALL iom_put ('snowvol', zcmo(:,:,40) ) ! field 40: snow volume 348 CALL iom_put ('sfxmec', zcmo(:,:,41) * rday ) ! field 41: salt flux from ridging rafting 349 CALL iom_put ('sfxres', zcmo(:,:,42) * rday ) ! field 42: salt flux from limupdate (resultant) 350 CALL iom_put ('icetrp', zcmo(:,:,43) ) ! field 43: ice volume transport 287 351 288 352 !----------------------------- … … 303 367 DO jj = 1, jpj 304 368 DO ji = 1, jpi 305 zinda = MAX( zzero , SIGN( zone , a_i(ji,jj,jl) - 1.0e-6 ) )306 zoi(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , 1.0e-6 ) * zinda369 zinda = MAX( zzero , SIGN( zone , a_i(ji,jj,jl) - epsi06 ) ) 370 zoi(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi06 ) * zinda 307 371 END DO 308 372 END DO … … 315 379 DO jj = 1, jpj 316 380 DO ji = 1, jpi 317 zinda = MAX( zzero , SIGN( zone , a_i(ji,jj,jl) - 1.0e-6 ) )381 zinda = MAX( zzero , SIGN( zone , a_i(ji,jj,jl) - epsi06 ) ) 318 382 zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0* & 319 ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rtt ), - 1.0e-6 ) ) * &383 ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rtt ), - epsi06 ) ) * & 320 384 zinda / nlay_i 321 385 END DO … … 349 413 CALL wrk_dealloc( jpi, jpj, jpnoumax, zcmo, zcmoa ) 350 414 CALL wrk_dealloc( jpi, jpj, jpl, zmaskitd, zoi, zei ) 415 416 IF( nn_timing == 1 ) CALL timing_stop('limwri') 351 417 352 418 END SUBROUTINE lim_wri … … 382 448 field_25, field_26, field_27, field_28, field_29, field_30, & 383 449 field_31, field_32, field_33, field_34, field_35, field_36, & 384 field_37 450 field_37, field_38, field_39, field_40, field_41, field_42, field_43 385 451 386 452 TYPE(FIELD) , DIMENSION(jpnoumax) :: zfield … … 393 459 field_25, field_26, field_27, field_28, field_29, field_30, & 394 460 field_31, field_32, field_33, field_34, field_35, field_36, & 395 field_37, add_diag_swi461 field_37, field_38, field_39, field_40, field_41, field_42, field_43, add_diag_swi 396 462 !!------------------------------------------------------------------- 397 463 … … 436 502 zfield(36) = field_36 437 503 zfield(37) = field_37 504 zfield(38) = field_38 505 zfield(39) = field_39 506 zfield(40) = field_40 507 zfield(41) = field_41 508 zfield(42) = field_42 509 zfield(43) = field_43 438 510 439 511 DO nf = 1, noumef … … 461 533 ! 462 534 END SUBROUTINE lim_wri_init 535 536 SUBROUTINE lim_wri_state( kt, kid, kh_i ) 537 !!--------------------------------------------------------------------- 538 !! *** ROUTINE lim_wri_state *** 539 !! 540 !! ** Purpose : create a NetCDF file named cdfile_name which contains 541 !! the instantaneous ice state and forcing fields for ice model 542 !! Used to find errors in the initial state or save the last 543 !! ocean state in case of abnormal end of a simulation 544 !! 545 !! History : 546 !! 4.1 ! 2013-06 (C. Rousset) 547 !!---------------------------------------------------------------------- 548 INTEGER, INTENT( in ) :: kt ! ocean time-step index) 549 INTEGER, INTENT( in ) :: kid , kh_i 550 !!---------------------------------------------------------------------- 551 !CALL histvert( kid, "icethi", "L levels","m", jpl , hi_mean , nz ) 552 553 CALL histdef( kid, "iicethic", "Ice thickness" , "m" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 554 CALL histdef( kid, "iiceconc", "Ice concentration" , "%" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 555 CALL histdef( kid, "iicetemp", "Ice temperature" , "C" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 556 CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 557 CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 558 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 ) 559 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 ) 560 CALL histdef( kid, "iicesflx", "Solar flux over ocean" , "w/m2" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 561 CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 562 CALL histdef( kid, "isnowpre", "Snow precipitation" , "kg/m2/s", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 563 CALL histdef( kid, "iicesali", "Ice salinity" , "PSU" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 564 CALL histdef( kid, "iicevolu", "Ice volume" , "m" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 565 CALL histdef( kid, "iicedive", "Ice divergence" , "10-8s-1", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 566 567 !CALL histdef( kid, "iice_itd", "Ice concentration by cat", "%" , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 568 !CALL histdef( kid, "iice_hid", "Ice thickness by cat" , "m" , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 569 !CALL histdef( kid, "iice_hsd", "Snow thickness by cat" , "m" , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 570 !CALL histdef( kid, "iice_std", "Ice salinity by cat" , "PSU" , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 571 572 CALL histend( kid, snc4set ) ! end of the file definition 573 574 CALL histwrite( kid, "iicethic", kt, icethi , jpi*jpj, (/1/) ) 575 CALL histwrite( kid, "iiceconc", kt, at_i , jpi*jpj, (/1/) ) 576 CALL histwrite( kid, "iicetemp", kt, tm_i - rtt , jpi*jpj, (/1/) ) 577 CALL histwrite( kid, "iicevelu", kt, u_ice , jpi*jpj, (/1/) ) 578 CALL histwrite( kid, "iicevelv", kt, v_ice , jpi*jpj, (/1/) ) 579 CALL histwrite( kid, "iicestru", kt, utau_ice , jpi*jpj, (/1/) ) 580 CALL histwrite( kid, "iicestrv", kt, vtau_ice , jpi*jpj, (/1/) ) 581 CALL histwrite( kid, "iicesflx", kt, qsr , jpi*jpj, (/1/) ) 582 CALL histwrite( kid, "iicenflx", kt, qns , jpi*jpj, (/1/) ) 583 CALL histwrite( kid, "isnowpre", kt, sprecip , jpi*jpj, (/1/) ) 584 CALL histwrite( kid, "iicesali", kt, smt_i , jpi*jpj, (/1/) ) 585 CALL histwrite( kid, "iicevolu", kt, vt_i , jpi*jpj, (/1/) ) 586 CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8 , jpi*jpj, (/1/) ) 587 588 !CALL histwrite( kid, "iice_itd", kt, a_i , jpi*jpj*jpl, (/1/) ) ! area 589 !CALL histwrite( kid, "iice_hid", kt, ht_i , jpi*jpj*jpl, (/1/) ) ! thickness 590 !CALL histwrite( kid, "iice_hsd", kt, ht_s , jpi*jpj*jpl, (/1/) ) ! snow depth 591 !CALL histwrite( kid, "iice_std", kt, sm_i , jpi*jpj*jpl, (/1/) ) ! salinity 592 593 END SUBROUTINE lim_wri_state 463 594 464 595 #else -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r3625 r4045 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) … … 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 … … 138 140 139 141 !!---------------------------------------------------------------------- 140 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)142 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 141 143 !! $Id$ 142 144 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 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_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
r3625 r4045 5 5 !!====================================================================== 6 6 !! History : 3.3 ! 2010-09 (M. Leclair) Original code 7 !! ! 2012-10 (C. Rousset) add iom_put 7 8 !!---------------------------------------------------------------------- 8 9 … … 21 22 USE bdy_par ! (for lk_bdy) 22 23 USE timing ! preformance summary 24 USE iom ! I/O manager 25 USE lib_fortran ! glob_sum 26 USE restart ! ocean restart 27 USE wrk_nemo ! work arrays 23 28 24 29 IMPLICIT NONE … … 26 31 27 32 PUBLIC dia_hsb ! routine called by step.F90 28 PUBLIC dia_hsb_init ! routine called by opa.F90 33 PUBLIC dia_hsb_init ! routine called by nemogcm.F90 34 PUBLIC dia_hsb_rst ! routine called by step.F90 29 35 30 36 LOGICAL, PUBLIC :: ln_diahsb = .FALSE. !: check the heat and salt budgets 31 37 32 INTEGER :: numhsb ! 33 REAL(dp) :: surf_tot , vol_tot ! 34 REAL(dp) :: frc_t , frc_s , frc_v ! global forcing trends 35 REAL(dp) :: fact1 ! conversion factors 36 REAL(dp) :: fact21 , fact22 ! - - 37 REAL(dp) :: fact31 , fact32 ! - - 38 REAL(dp), DIMENSION(:,:) , ALLOCATABLE :: surf , ssh_ini ! 39 REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: hc_loc_ini, sc_loc_ini, e3t_ini ! 38 REAL(wp), SAVE :: frc_t , frc_s , frc_v ! global forcing trends 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssh_ini ! 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hc_loc_ini, sc_loc_ini, e3t_ini ! 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hcssh_loc_ini, scssh_loc_ini ! 40 42 41 43 !! * Substitutions … … 61 63 !! - Compute the contribution of forcing and remove it from these deviations 62 64 !! 63 !! ** Action : Write the results in the 'heat_salt_volume_budgets.txt' ASCII file64 65 !!--------------------------------------------------------------------------- 65 66 INTEGER, INTENT(in) :: kt ! ocean time-step index 66 67 !! 67 68 INTEGER :: jk ! dummy loop indice 68 REAL(dp) :: zdiff_hc , zdiff_sc ! heat and salt content variations 69 REAL(dp) :: zdiff_v1 , zdiff_v2 ! volume variation 70 REAL(dp) :: z1_rau0 ! local scalars 71 REAL(dp) :: zdeltat ! - - 72 REAL(dp) :: z_frc_trd_t , z_frc_trd_s ! - - 73 REAL(dp) :: z_frc_trd_v ! - - 74 !!--------------------------------------------------------------------------- 75 IF( nn_timing == 1 ) CALL timing_start('dia_hsb') 76 69 REAL(wp) :: zdiff_hc , zdiff_sc ! heat and salt content variations 70 REAL(wp) :: zdiff_v1 , zdiff_v2 ! volume variation 71 REAL(wp) :: z_hc , z_sc ! heat and salt content 72 REAL(wp) :: z_v1 , z_v2 ! volume 73 REAL(wp) :: z1_rau0 ! local scalars 74 REAL(wp) :: zdeltat ! - - 75 REAL(wp) :: z_frc_trd_t , z_frc_trd_s ! - - 76 REAL(wp) :: z_frc_trd_v ! - - 77 REAL(wp), POINTER, DIMENSION(:,:) :: zsurf ! 78 !!--------------------------------------------------------------------------- 79 IF( nn_timing == 1 ) CALL timing_start('dia_hsb') 80 81 CALL wrk_alloc( jpi, jpj, zsurf ) 82 83 zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:) ! masked surface grid cell area 84 77 85 ! ------------------------- ! 78 86 ! 1 - Trends due to forcing ! 79 87 ! ------------------------- ! 80 88 z1_rau0 = 1.e0 / rau0 81 z_frc_trd_v = z1_rau0 * SUM( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) )! volume fluxes82 z_frc_trd_t = SUM( sbc_tsc(:,:,jp_tem) * surf(:,:) )! heat fluxes83 z_frc_trd_s = SUM( sbc_tsc(:,:,jp_sal) * surf(:,:) )! salt fluxes89 z_frc_trd_v = z1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * zsurf(:,:) ) ! volume fluxes 90 z_frc_trd_t = glob_sum( sbc_tsc(:,:,jp_tem) * zsurf(:,:) ) ! heat fluxes 91 z_frc_trd_s = glob_sum( sbc_tsc(:,:,jp_sal) * zsurf(:,:) ) ! salt fluxes 84 92 ! Add penetrative solar radiation 85 IF( ln_traqsr ) z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * SUM( qsr (:,:) *surf(:,:) )93 IF( ln_traqsr ) z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr (:,:) * zsurf(:,:) ) 86 94 ! Add geothermal heat flux 87 IF( ln_trabbc ) z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * SUM( qgh_trd0(:,:) * surf(:,:) ) 88 IF( lk_mpp ) THEN 89 CALL mpp_sum( z_frc_trd_v ) 90 CALL mpp_sum( z_frc_trd_t ) 91 ENDIF 95 IF( ln_trabbc ) z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qgh_trd0(:,:) * zsurf(:,:) ) 96 ! 92 97 frc_v = frc_v + z_frc_trd_v * rdt 93 98 frc_t = frc_t + z_frc_trd_t * rdt 94 99 frc_s = frc_s + z_frc_trd_s * rdt 95 100 96 ! ----------------------- !97 ! 2 - Content variations !98 ! ----------------------- !99 zdiff_v2 = 0. d0100 zdiff_hc = 0. d0101 zdiff_sc = 0. d0101 ! ------------------------ ! 102 ! 2a - Content variations ! 103 ! ------------------------ ! 104 zdiff_v2 = 0._wp 105 zdiff_hc = 0._wp 106 zdiff_sc = 0._wp 102 107 ! volume variation (calculated with ssh) 103 zdiff_v1 = SUM( surf(:,:) * tmask(:,:,1) * ( sshn(:,:) - ssh_ini(:,:) ) )108 zdiff_v1 = glob_sum( zsurf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) ) 104 109 DO jk = 1, jpkm1 105 110 ! volume variation (calculated with scale factors) 106 zdiff_v2 = zdiff_v2 + SUM( surf(:,:) * tmask(:,:,jk) & 107 & * ( fse3t_n(:,:,jk) & 108 & - e3t_ini(:,:,jk) ) ) 111 zdiff_v2 = zdiff_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) ) 109 112 ! heat content variation 110 zdiff_hc = zdiff_hc + SUM( surf(:,:) * tmask(:,:,jk) & 111 & * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) & 113 zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) & 112 114 & - hc_loc_ini(:,:,jk) ) ) 113 115 ! salt content variation 114 zdiff_sc = zdiff_sc + SUM( surf(:,:) * tmask(:,:,jk) & 115 & * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) & 116 zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) & 116 117 & - sc_loc_ini(:,:,jk) ) ) 117 118 ENDDO 118 119 119 IF( lk_mpp ) THEN120 CALL mpp_sum( zdiff_hc )121 CALL mpp_sum( zdiff_sc )122 CALL mpp_sum( zdiff_v1 )123 CALL mpp_sum( zdiff_v2 )124 ENDIF125 126 120 ! Substract forcing from heat content, salt content and volume variations 127 zdiff_v1 = zdiff_v1 - frc_v 128 zdiff_v2 = zdiff_v2 - frc_v 129 zdiff_hc = zdiff_hc - frc_t 130 zdiff_sc = zdiff_sc - frc_s 121 !frc_v = zdiff_v2 - frc_v 122 !frc_t = zdiff_hc - frc_t 123 !frc_s = zdiff_sc - frc_s 131 124 125 ! add ssh if not vvl 126 #if ! defined key_vvl 127 zdiff_v2 = zdiff_v2 + zdiff_v1 128 zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_tem) & 129 & - hcssh_loc_ini(:,:) ) ) 130 zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_sal) & 131 & - scssh_loc_ini(:,:) ) ) 132 #endif 133 ! 134 ! ----------------------- ! 135 ! 2b - Content ! 136 ! ----------------------- ! 137 z_v2 = 0._wp 138 z_hc = 0._wp 139 z_sc = 0._wp 140 ! volume (calculated with ssh) 141 z_v1 = glob_sum( zsurf(:,:) * sshn(:,:) ) 142 DO jk = 1, jpkm1 143 ! volume (calculated with scale factors) 144 z_v2 = z_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) ) 145 ! heat content 146 z_hc = z_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) ) 147 ! salt content 148 z_sc = z_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) ) 149 ENDDO 150 ! add ssh if not vvl 151 #if ! defined key_vvl 152 z_v2 = z_v2 + z_v1 153 z_hc = z_hc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_tem) ) 154 z_sc = z_sc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_sal) ) 155 #endif 156 132 157 ! ----------------------- ! 133 158 ! 3 - Diagnostics writing ! 134 159 ! ----------------------- ! 135 160 zdeltat = 1.e0 / ( ( kt - nit000 + 1 ) * rdt ) 136 WRITE(numhsb , 9020) kt , zdiff_hc / vol_tot , zdiff_hc * fact1 * zdeltat, & 137 & zdiff_sc / vol_tot , zdiff_sc * fact21 * zdeltat, zdiff_sc * fact22 * zdeltat, & 138 & zdiff_v1 , zdiff_v1 * fact31 * zdeltat, zdiff_v1 * fact32 * zdeltat, & 139 & zdiff_v2 , zdiff_v2 * fact31 * zdeltat, zdiff_v2 * fact32 * zdeltat 140 141 IF ( kt == nitend ) CLOSE( numhsb ) 142 161 ! 162 CALL iom_put( 'bgtemper',z_hc / z_v2 ) ! Temperature (C) 163 CALL iom_put( 'bgsaline',z_sc / z_v2 ) ! Salinity (psu) 164 !CALL iom_put( 'bgheatco',zdiff_hc*fact1*zdeltat ) ! Equivalent heat flux (W/m2) 165 !CALL iom_put( 'bgsaltco',zdiff_sc*fact21*zdeltat ) ! Equivalent water flux (mm/s) 166 CALL iom_put( 'bgheatco',zdiff_hc * rau0 * rcp * 1.e-9_wp ) ! Heat content variation (10^9 J) 167 CALL iom_put( 'bgsaltco',zdiff_sc * 1.e-9 ) ! Salt content variation (psu*km3) 168 CALL iom_put( 'bgvolssh',zdiff_v1 * 1.e-9 ) ! volume ssh (km3) 169 CALL iom_put( 'bgsshtot',zdiff_v1 / glob_sum(zsurf) ) ! ssh (m) 170 CALL iom_put( 'bgvoltot',zdiff_v2 * 1.e-9 ) ! volume total (km3) 171 CALL iom_put( 'bgfrcvol',frc_v * 1.e-9 ) ! vol - surface forcing (volume) 172 CALL iom_put( 'bgfrctem',frc_t * rau0 * rcp * 1.e-9_wp ) ! hc - surface forcing (heat content) 173 CALL iom_put( 'bgfrcsal',frc_s * 1.e-9 ) ! sc - surface forcing (salt content) 174 ! 175 CALL wrk_dealloc( jpi, jpj, zsurf ) 176 ! 143 177 IF( nn_timing == 1 ) CALL timing_stop('dia_hsb') 144 145 9020 FORMAT(I5,11D15.7) 146 ! 178 ! 147 179 END SUBROUTINE dia_hsb 148 180 … … 160 192 !! - Compute coefficients for conversion 161 193 !!--------------------------------------------------------------------------- 162 CHARACTER (len=32) :: cl_name ! output file name163 194 INTEGER :: jk ! dummy loop indice 164 195 INTEGER :: ierror ! local integer … … 180 211 IF( .NOT. ln_diahsb ) RETURN 181 212 182 ! ------------------- ! 183 ! 1 - Allocate memory ! 184 ! ------------------- ! 185 ALLOCATE( hc_loc_ini(jpi,jpj,jpk), STAT=ierror ) 186 IF( ierror > 0 ) THEN 187 CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) ; RETURN 188 ENDIF 189 ALLOCATE( sc_loc_ini(jpi,jpj,jpk), STAT=ierror ) 190 IF( ierror > 0 ) THEN 191 CALL ctl_stop( 'dia_hsb: unable to allocate sc_loc_ini' ) ; RETURN 192 ENDIF 193 ALLOCATE( e3t_ini(jpi,jpj,jpk) , STAT=ierror ) 194 IF( ierror > 0 ) THEN 195 CALL ctl_stop( 'dia_hsb: unable to allocate e3t_ini' ) ; RETURN 196 ENDIF 197 ALLOCATE( surf(jpi,jpj) , STAT=ierror ) 198 IF( ierror > 0 ) THEN 199 CALL ctl_stop( 'dia_hsb: unable to allocate surf' ) ; RETURN 200 ENDIF 201 ALLOCATE( ssh_ini(jpi,jpj) , STAT=ierror ) 202 IF( ierror > 0 ) THEN 203 CALL ctl_stop( 'dia_hsb: unable to allocate ssh_ini' ) ; RETURN 204 ENDIF 205 206 ! ----------------------------------------------- ! 207 ! 2 - Time independant variables and file opening ! 208 ! ----------------------------------------------- ! 209 WRITE(numout,*) "dia_hsb: heat salt volume budgets activated" 210 WRITE(numout,*) "~~~~~~~ output written in the 'heat_salt_volume_budgets.txt' ASCII file" 211 IF( lk_obc .or. lk_bdy ) THEN 212 CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' ) 213 ENDIF 214 cl_name = 'heat_salt_volume_budgets.txt' ! name of output file 215 surf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:) ! masked surface grid cell area 216 surf_tot = SUM( surf(:,:) ) ! total ocean surface area 217 vol_tot = 0.d0 ! total ocean volume 218 DO jk = 1, jpkm1 219 vol_tot = vol_tot + SUM( surf(:,:) * tmask(:,:,jk) & 220 & * fse3t_n(:,:,jk) ) 221 END DO 222 IF( lk_mpp ) THEN 223 CALL mpp_sum( vol_tot ) 224 CALL mpp_sum( surf_tot ) 225 ENDIF 226 227 CALL ctl_opn( numhsb , cl_name , 'UNKNOWN' , 'FORMATTED' , 'SEQUENTIAL' , 1 , numout , lwp , 1 ) 228 ! 12345678901234567890123456789012345678901234567890123456789012345678901234567890 -> 80 229 WRITE( numhsb, 9010 ) "kt | heat content budget | salt content budget ", & 230 ! 123456789012345678901234567890123456789012345 -> 45 231 & "| volume budget (ssh) ", & 232 ! 678901234567890123456789012345678901234567890 -> 45 233 & "| volume budget (e3t) " 234 WRITE( numhsb, 9010 ) " | [C] [W/m2] | [psu] [mmm/s] [SV] ", & 235 & "| [m3] [mmm/s] [SV] ", & 236 & "| [m3] [mmm/s] [SV] " 237 238 ! --------------- ! 239 ! 3 - Conversions ! (factors will be multiplied by duration afterwards) 240 ! --------------- ! 241 242 ! heat content variation => equivalent heat flux: 243 fact1 = rau0 * rcp / surf_tot ! [C*m3] -> [W/m2] 244 ! salt content variation => equivalent EMP and equivalent "flow": 245 fact21 = 1.e3 / ( soce * surf_tot ) ! [psu*m3] -> [mm/s] 246 fact22 = 1.e-6 / soce ! [psu*m3] -> [Sv] 247 ! volume variation => equivalent EMP and equivalent "flow": 248 fact31 = 1.e3 / surf_tot ! [m3] -> [mm/s] 249 fact32 = 1.e-6 ! [m3] -> [SV] 250 251 ! ---------------------------------- ! 252 ! 4 - initial conservation variables ! 253 ! ---------------------------------- ! 254 ssh_ini(:,:) = sshn(:,:) ! initial ssh 255 DO jk = 1, jpk 256 e3t_ini (:,:,jk) = fse3t_n(:,:,jk) ! initial vertical scale factors 257 hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) ! initial heat content 258 sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) ! initial salt content 259 END DO 260 frc_v = 0.d0 ! volume trend due to forcing 261 frc_t = 0.d0 ! heat content - - - - 262 frc_s = 0.d0 ! salt content - - - - 263 ! 264 9010 FORMAT(A80,A45,A45) 213 ! ------------------- ! 214 ! 1 - Allocate memory ! 215 ! ------------------- ! 216 ALLOCATE( hc_loc_ini(jpi,jpj,jpk), STAT=ierror ) 217 IF( ierror > 0 ) THEN 218 CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) ; RETURN 219 ENDIF 220 ALLOCATE( sc_loc_ini(jpi,jpj,jpk), STAT=ierror ) 221 IF( ierror > 0 ) THEN 222 CALL ctl_stop( 'dia_hsb: unable to allocate sc_loc_ini' ) ; RETURN 223 ENDIF 224 ALLOCATE( hcssh_loc_ini(jpi,jpj), STAT=ierror ) 225 IF( ierror > 0 ) THEN 226 CALL ctl_stop( 'dia_hsb: unable to allocate hcssh_loc_ini' ) ; RETURN 227 ENDIF 228 ALLOCATE( scssh_loc_ini(jpi,jpj), STAT=ierror ) 229 IF( ierror > 0 ) THEN 230 CALL ctl_stop( 'dia_hsb: unable to allocate scssh_loc_ini' ) ; RETURN 231 ENDIF 232 ALLOCATE( e3t_ini(jpi,jpj,jpk) , STAT=ierror ) 233 IF( ierror > 0 ) THEN 234 CALL ctl_stop( 'dia_hsb: unable to allocate e3t_ini' ) ; RETURN 235 ENDIF 236 ALLOCATE( ssh_ini(jpi,jpj) , STAT=ierror ) 237 IF( ierror > 0 ) THEN 238 CALL ctl_stop( 'dia_hsb: unable to allocate ssh_ini' ) ; RETURN 239 ENDIF 240 241 ! ----------------------------------------------- ! 242 ! 2 - Time independant variables and file opening ! 243 ! ----------------------------------------------- ! 244 WRITE(numout,*) "dia_hsb: heat salt volume budgets activated" 245 IF( lk_obc .or. lk_bdy ) THEN 246 CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' ) 247 ENDIF 248 249 ! ---------------------------------- ! 250 ! 4 - initial conservation variables ! 251 ! ---------------------------------- ! 252 !ssh_ini(:,:) = sshn(:,:) ! initial ssh 253 !DO jk = 1, jpk 254 ! e3t_ini (:,:,jk) = fse3t_n(:,:,jk) ! initial vertical scale factors 255 ! hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) ! initial heat content 256 ! sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) ! initial salt content 257 !END DO 258 !hcssh_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) ! initial heat content in ssh 259 !scssh_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:) ! initial salt content in ssh 260 !frc_v = 0._wp ! volume trend due to forcing 261 !frc_t = 0._wp ! heat content - - - - 262 !frc_s = 0._wp ! salt content - - - - 263 ! 264 CALL dia_hsb_rst( nit000, 'READ' ) !* read or initialize all required files 265 265 ! 266 266 END SUBROUTINE dia_hsb_init 267 268 SUBROUTINE dia_hsb_rst( kt, cdrw ) 269 !!--------------------------------------------------------------------- 270 !! *** ROUTINE limdia_rst *** 271 !! 272 !! ** Purpose : Read or write DIA file in restart file 273 !! 274 !! ** Method : use of IOM library 275 !!---------------------------------------------------------------------- 276 INTEGER , INTENT(in) :: kt ! ocean time-step 277 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 278 ! 279 INTEGER :: jk ! 280 INTEGER :: id1 ! local integers 281 !!---------------------------------------------------------------------- 282 ! 283 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 284 IF( ln_rstart ) THEN !* Read the restart file 285 !id1 = iom_varid( numror, 'frc_vol' , ldstop = .FALSE. ) 286 ! 287 CALL iom_get( numror, 'frc_v', frc_v ) 288 CALL iom_get( numror, 'frc_t', frc_t ) 289 CALL iom_get( numror, 'frc_s', frc_s ) 290 291 CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini ) 292 CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini ) 293 CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini ) 294 CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini ) 295 CALL iom_get( numror, jpdom_autoglo, 'hcssh_loc_ini', hcssh_loc_ini ) 296 CALL iom_get( numror, jpdom_autoglo, 'scssh_loc_ini', scssh_loc_ini ) 297 ELSE 298 ssh_ini(:,:) = sshn(:,:) ! initial ssh 299 DO jk = 1, jpk 300 e3t_ini (:,:,jk) = fse3t_n(:,:,jk) ! initial vertical scale factors 301 hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) ! initial heat content 302 sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) ! initial salt content 303 END DO 304 hcssh_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) ! initial heat content in ssh 305 scssh_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:) ! initial salt content in ssh 306 frc_v = 0._wp 307 frc_t = 0._wp 308 frc_s = 0._wp 309 ENDIF 310 311 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 312 ! ! ------------------- 313 IF(lwp) WRITE(numout,*) '---- dia-rst ----' 314 CALL iom_rstput( kt, nitrst, numrow, 'frc_v' , frc_v ) 315 CALL iom_rstput( kt, nitrst, numrow, 'frc_t' , frc_t ) 316 CALL iom_rstput( kt, nitrst, numrow, 'frc_s' , frc_s ) 317 318 CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini ) 319 CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini ) 320 CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini ) 321 CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini ) 322 CALL iom_rstput( kt, nitrst, numrow, 'hcssh_loc_ini', hcssh_loc_ini ) 323 CALL iom_rstput( kt, nitrst, numrow, 'scssh_loc_ini', scssh_loc_ini ) 324 ! 325 ENDIF 326 ! 327 END SUBROUTINE dia_hsb_rst 267 328 268 329 !!====================================================================== -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r3704 r4045 48 48 #if defined key_lim2 49 49 USE limwri_2 50 #elif defined key_lim3 51 USE limwri 50 52 #endif 51 53 USE lib_mpp ! MPP library … … 842 844 #if defined key_lim2 843 845 CALL lim_wri_state_2( kt, id_i, nh_i ) 846 #elif defined key_lim3 847 CALL lim_wri_state( kt, id_i, nh_i ) 844 848 #else 845 849 CALL histend( id_i, snc4chunks=snc4set ) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
r3983 r4045 392 392 393 393 394 FUNCTION iom_varid ( kiomid, cdvar, kdimsz, ldstop )394 FUNCTION iom_varid ( kiomid, cdvar, kdimsz, kndims, ldstop ) 395 395 !!----------------------------------------------------------------------- 396 396 !! *** FUNCTION iom_varid *** … … 401 401 CHARACTER(len=*) , INTENT(in ) :: cdvar ! name of the variable 402 402 INTEGER, DIMENSION(:), INTENT( out), OPTIONAL :: kdimsz ! size of the dimensions 403 INTEGER, INTENT( out), OPTIONAL :: kndims ! size of the dimensions 403 404 LOGICAL , INTENT(in ), OPTIONAL :: ldstop ! stop if looking for non-existing variable (default = .TRUE.) 404 405 ! … … 432 433 SELECT CASE (iom_file(kiomid)%iolib) 433 434 CASE (jpioipsl ) ; iom_varid = iom_ioipsl_varid( kiomid, cdvar, iiv, kdimsz ) 434 CASE (jpnf90 ) ; iom_varid = iom_nf90_varid ( kiomid, cdvar, iiv, kdimsz )435 CASE (jpnf90 ) ; iom_varid = iom_nf90_varid ( kiomid, cdvar, iiv, kdimsz, kndims ) 435 436 CASE (jprstdimg) ; iom_varid = -1 ! all variables are listed in iom_file 436 437 CASE DEFAULT … … 453 454 ENDIF 454 455 ENDIF 456 IF( PRESENT(kndims) ) kndims = iom_file(kiomid)%ndims(iiv) 455 457 ENDIF 456 458 ENDIF … … 1185 1187 WRITE(cl1,'(i1)') 1 ; CALL iom_set_field_attr('field_definition', freq_op = cl1//'ts', freq_offset='0ts') 1186 1188 WRITE(cl1,'(i1)') nn_fsbc ; CALL iom_set_field_attr('SBC' , freq_op = cl1//'ts', freq_offset='0ts') 1189 WRITE(cl1,'(i1)') nn_fsbc ; CALL iom_set_field_attr('SBC_scalar' , freq_op = cl1//'ts', freq_offset='0ts') 1187 1190 WRITE(cl1,'(i1)') nn_dttrc ; CALL iom_set_field_attr('ptrc_T' , freq_op = cl1//'ts', freq_offset='0ts') 1188 1191 WRITE(cl1,'(i1)') nn_dttrc ; CALL iom_set_field_attr('diad_T' , freq_op = cl1//'ts', freq_offset='0ts') -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90
r2715 r4045 181 181 182 182 183 FUNCTION iom_nf90_varid ( kiomid, cdvar, kiv, kdimsz )183 FUNCTION iom_nf90_varid ( kiomid, cdvar, kiv, kdimsz, kndims ) 184 184 !!----------------------------------------------------------------------- 185 185 !! *** FUNCTION iom_varid *** … … 191 191 INTEGER , INTENT(in ) :: kiv ! 192 192 INTEGER, DIMENSION(:), INTENT( out), OPTIONAL :: kdimsz ! size of the dimensions 193 INTEGER, INTENT( out), OPTIONAL :: kndims ! size of the dimensions 193 194 ! 194 195 INTEGER :: iom_nf90_varid ! iom variable Id … … 242 243 ENDIF 243 244 ENDIF 245 IF( PRESENT(kndims) ) kndims = iom_file(kiomid)%ndims(kiv) 244 246 ELSE 245 247 iom_nf90_varid = -1 ! variable not found, return error code: -1
Note: See TracChangeset
for help on using the changeset viewer.