Changeset 4972 for branches/2014/dev_MERGE_2014
- Timestamp:
- 2014-12-05T16:28:42+01:00 (9 years ago)
- Location:
- branches/2014/dev_MERGE_2014/NEMOGCM
- Files:
-
- 23 edited
-
CONFIG/SHARED/namelist_ice_lim3_ref (modified) (1 diff)
-
CONFIG/cfg.txt (modified) (1 diff)
-
NEMO/LIM_SRC_3/ice.F90 (modified) (9 diffs)
-
NEMO/LIM_SRC_3/limadv.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limdiahsb.F90 (modified) (6 diffs)
-
NEMO/LIM_SRC_3/limdyn.F90 (modified) (11 diffs)
-
NEMO/LIM_SRC_3/limitd_me.F90 (modified) (2 diffs)
-
NEMO/LIM_SRC_3/limitd_th.F90 (modified) (6 diffs)
-
NEMO/LIM_SRC_3/limrhg.F90 (modified) (36 diffs)
-
NEMO/LIM_SRC_3/limrst.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limsbc.F90 (modified) (3 diffs)
-
NEMO/LIM_SRC_3/limthd.F90 (modified) (12 diffs)
-
NEMO/LIM_SRC_3/limthd_dh.F90 (modified) (15 diffs)
-
NEMO/LIM_SRC_3/limthd_dif.F90 (modified) (25 diffs)
-
NEMO/LIM_SRC_3/limthd_ent.F90 (modified) (3 diffs)
-
NEMO/LIM_SRC_3/limthd_lac.F90 (modified) (9 diffs)
-
NEMO/LIM_SRC_3/limtrp.F90 (modified) (9 diffs)
-
NEMO/LIM_SRC_3/limupdate1.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limupdate2.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limvar.F90 (modified) (19 diffs)
-
NEMO/LIM_SRC_3/limwri.F90 (modified) (13 diffs)
-
NEMO/OPA_SRC/BDY/bdyice_lim.F90 (modified) (11 diffs)
-
NEMO/OPA_SRC/SBC/sbcice_lim.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_MERGE_2014/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r4924 r4972 46 46 om = 0.5 ! relaxation constant 47 47 cw = 5.0e-03 ! drag coefficient for oceanic stress 48 angvg = 0.0 ! turning angle for oceanic stress49 48 pstar = 2.0e+04 ! 1st bulk-rheology parameter 50 49 c_rhg = 20.0 ! 2nd bulk-rhelogy parameter -
branches/2014/dev_MERGE_2014/NEMOGCM/CONFIG/cfg.txt
r4946 r4972 4 4 ORCA2_SAS_LIM OPA_SRC SAS_SRC LIM_SRC_2 NST_SRC 5 5 C1D_PAPA OPA_SRC 6 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC7 6 GYRE_BFM OPA_SRC TOP_SRC 8 7 AMM12 OPA_SRC 9 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC10 8 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 11 9 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 12 10 ISOMIP OPA_SRC 13 11 GYRE OPA_SRC 12 SPITZ025 OPA_SRC LIM_SRC_3 13 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 14 GYRE_LONG OPA_SRC 15 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 16 GYRE_4 OPA_SRC 17 ORCA2LIMPIS_LONG OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r4924 r4972 170 170 REAL(wp), PUBLIC :: om !: relaxation constant 171 171 REAL(wp), PUBLIC :: cw !: drag coefficient for oceanic stress 172 REAL(wp), PUBLIC :: angvg !: turning angle for oceanic stress173 172 REAL(wp), PUBLIC :: pstar !: determines ice strength (N/M), Hibler JPO79 174 173 REAL(wp), PUBLIC :: c_rhg !: determines changes in ice strength … … 179 178 REAL(wp), PUBLIC :: relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 180 179 REAL(wp), PUBLIC :: alphaevp !: coeficient of the internal stresses 181 REAL(wp), PUBLIC :: unit_fac = 1.e+09_wp !: conversion factor for ice / snow enthalpy182 180 REAL(wp), PUBLIC :: hminrhg !: ice volume (a*h, in m) below which ice velocity is set to ocean velocity 183 181 … … 224 222 REAL(wp), PUBLIC :: usecc2 !: = 1.0 / ( ecc * ecc ) 225 223 REAL(wp), PUBLIC :: rhoco !: = rau0 * cw 226 REAL(wp), PUBLIC :: sangvg, cangvg !: sin and cos of the turning angle for ocean stress 227 REAL(wp), PUBLIC :: pstarh !: pstar / 2.0 224 225 ! !!** switch for presence of ice or not 226 REAL(wp), PUBLIC :: rswitch 227 228 ! !!** define some parameters 229 REAL(wp), PUBLIC, PARAMETER :: unit_fac = 1.e+09_wp !: conversion factor for ice / snow enthalpy 230 REAL(wp), PUBLIC, PARAMETER :: epsi06 = 1.e-06_wp !: small number 231 REAL(wp), PUBLIC, PARAMETER :: epsi10 = 1.e-10_wp !: small number 232 REAL(wp), PUBLIC, PARAMETER :: epsi20 = 1.e-20_wp !: small number 228 233 229 234 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_oce, v_oce !: surface ocean velocity used in ice dynamics … … 364 369 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_sm_i_fl , d_sm_i_gd !: 365 370 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_sm_i_se , d_sm_i_si , d_sm_i_la !: 366 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_oa_i_thd , d_oa_i_trp , s_i_newice!:371 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_oa_i_thd , d_oa_i_trp !: 367 372 368 373 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: d_e_s_thd , d_e_s_trp !: … … 395 400 LOGICAL , PUBLIC :: ln_limdiaout !: flag for ice diag (T) or not (F) 396 401 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dv_dt_thd !: thermodynamic growth rates 397 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: izero398 402 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vi !: transport of ice volume 399 403 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vs !: transport of snw volume … … 418 422 INTEGER :: ice_alloc 419 423 ! 420 INTEGER :: ierr( 20), ii424 INTEGER :: ierr(19), ii 421 425 !!----------------------------------------------------------------- 422 426 … … 448 452 & hfx_thd(jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) , STAT=ierr(ii) ) 449 453 450 ii = ii + 1451 ALLOCATE( dh_i_surf2D(jpi,jpj) , dh_i_bott2D(jpi,jpj) , q_s(jpi,jpj) , STAT=ierr(ii) )452 453 454 ! * Ice global state variables 454 455 ii = ii + 1 … … 497 498 & d_v_i_thd(jpi,jpj,jpl) , d_v_i_trp (jpi,jpj,jpl) , d_smv_i_thd(jpi,jpj,jpl) , d_smv_i_trp(jpi,jpj,jpl) , & 498 499 & d_sm_i_fl(jpi,jpj,jpl) , d_sm_i_gd (jpi,jpj,jpl) , d_sm_i_se (jpi,jpj,jpl) , d_sm_i_si (jpi,jpj,jpl) , & 499 & d_sm_i_la(jpi,jpj,jpl) , d_oa_i_thd(jpi,jpj,jpl) , d_oa_i_trp (jpi,jpj,jpl) , s_i_newice (jpi,jpj,jpl) ,&500 & d_sm_i_la(jpi,jpj,jpl) , d_oa_i_thd(jpi,jpj,jpl) , d_oa_i_trp (jpi,jpj,jpl) , & 500 501 & STAT=ierr(ii) ) 501 502 ii = ii + 1 502 ALLOCATE( d_e_s_thd(jpi,jpj,nlay_s,jpl) , d_e_i_thd(jpi,jpj,nlay_i +1,jpl) , d_u_ice_dyn(jpi,jpj) , &503 & d_e_s_trp(jpi,jpj,nlay_s,jpl) , d_e_i_trp(jpi,jpj,nlay_i +1,jpl) , d_v_ice_dyn(jpi,jpj) , STAT=ierr(ii) )503 ALLOCATE( d_e_s_thd(jpi,jpj,nlay_s,jpl) , d_e_i_thd(jpi,jpj,nlay_i,jpl) , d_u_ice_dyn(jpi,jpj) , & 504 & d_e_s_trp(jpi,jpj,nlay_s,jpl) , d_e_i_trp(jpi,jpj,nlay_i,jpl) , d_v_ice_dyn(jpi,jpj) , STAT=ierr(ii) ) 504 505 505 506 ! * Ice thickness distribution variables … … 509 510 ! * Ice diagnostics 510 511 ii = ii + 1 511 ALLOCATE( dv_dt_thd(jpi,jpj,jpl), izero (jpi,jpj,jpl),&512 ALLOCATE( dv_dt_thd(jpi,jpj,jpl), & 512 513 & diag_trp_vi(jpi,jpj), diag_trp_vs (jpi,jpj), diag_trp_ei(jpi,jpj), & 513 514 & diag_trp_es(jpi,jpj), diag_heat_dhc(jpi,jpj), STAT=ierr(ii) ) -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
r4924 r4972 30 30 PUBLIC lim_adv_x ! called by lim_trp 31 31 PUBLIC lim_adv_y ! called by lim_trp 32 33 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values34 32 35 33 !! * Substitutions -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
r4967 r4972 37 37 real(wp) :: frc_sal, frc_vol ! global forcing trends 38 38 real(wp) :: bg_grme ! global ice growth+melt trends 39 REAL(wp) :: epsi06 = 1.e-6_wp ! small number40 39 41 40 !! * Substitutions … … 68 67 real(wp) :: z_frc_vol, z_frc_sal, z_bg_grme 69 68 real(wp) :: z1_area ! - - 70 real(wp) :: zinda, zindb71 69 REAL(wp) :: ztmp 72 70 !!--------------------------------------------------------------------------- … … 78 76 z1_area = 1._wp / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 ) 79 77 80 zinda= MAX( 0._wp , SIGN( 1._wp , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) )78 rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) ) 81 79 ! ----------------------- ! 82 80 ! 1 - Content variations ! … … 92 90 93 91 ! Volume 94 ztmp = zinda* z1_area * r1_rau0 * rday92 ztmp = rswitch * z1_area * r1_rau0 * rday 95 93 zbg_vfx = ztmp * glob_sum( emp(:,:) * area(:,:) * tms(:,:) ) 96 94 zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * area(:,:) * tms(:,:) ) … … 155 153 ! 3 - Diagnostics writing ! 156 154 ! ----------------------- ! 157 zindb= MAX( 0._wp , SIGN( 1._wp , zbg_ivo - epsi06 ) )155 rswitch = MAX( 0._wp , SIGN( 1._wp , zbg_ivo - epsi06 ) ) 158 156 ! 159 157 IF( iom_use('ibgvoltot') ) & … … 164 162 CALL iom_put( 'ibgarea' , zbg_are * 1.e-6 ) ! ice area (km2) 165 163 IF( iom_use('ibgsaline') ) & 166 CALL iom_put( 'ibgsaline' , zindb* zbg_sal / MAX( zbg_ivo, epsi06 ) ) ! ice saline (psu)164 CALL iom_put( 'ibgsaline' , rswitch * zbg_sal / MAX( zbg_ivo, epsi06 ) ) ! ice saline (psu) 167 165 IF( iom_use('ibgtemper') ) & 168 CALL iom_put( 'ibgtemper' , zindb* zbg_tem / MAX( zbg_ivo, epsi06 ) ) ! ice temper (C)166 CALL iom_put( 'ibgtemper' , rswitch * zbg_tem / MAX( zbg_ivo, epsi06 ) ) ! ice temper (C) 169 167 CALL iom_put( 'ibgheatco' , zbg_ihc ) ! ice heat content (1.e20 J) 170 168 CALL iom_put( 'sbgheatco' , zbg_shc ) ! snw heat content (1.e20 J) -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r4924 r4972 64 64 INTEGER :: i_j1, i_jpj ! Starting/ending j-indices for rheology 65 65 REAL(wp) :: zcoef ! local scalar 66 REAL(wp), POINTER, DIMENSION(:) :: z ind! i-averaged indicator of sea-ice66 REAL(wp), POINTER, DIMENSION(:) :: zswitch ! i-averaged indicator of sea-ice 67 67 REAL(wp), POINTER, DIMENSION(:) :: zmsk ! i-averaged of tmask 68 68 REAL(wp), POINTER, DIMENSION(:,:) :: zu_io, zv_io ! ice-ocean velocity … … 74 74 75 75 CALL wrk_alloc( jpi, jpj, zu_io, zv_io ) 76 CALL wrk_alloc( jpj, z ind, zmsk )76 CALL wrk_alloc( jpj, zswitch, zmsk ) 77 77 78 78 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) … … 100 100 ! 101 101 DO jj = 1, jpj 102 z ind(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line102 zswitch(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line 103 103 zmsk(jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line 104 104 END DO … … 110 110 i_j1 = njeq 111 111 i_jpj = jpj 112 DO WHILE ( i_j1 <= jpj .AND. z ind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )112 DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 113 113 i_j1 = i_j1 + 1 114 114 END DO … … 120 120 i_j1 = 1 121 121 i_jpj = njeq 122 DO WHILE ( i_jpj >= 1 .AND. z ind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )122 DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 123 123 i_jpj = i_jpj - 1 124 124 END DO … … 132 132 ! ! latitude strip 133 133 i_j1 = 1 134 DO WHILE ( i_j1 <= jpj .AND. z ind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )134 DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 135 135 i_j1 = i_j1 + 1 136 136 END DO … … 138 138 139 139 i_jpj = jpj 140 DO WHILE ( i_jpj >= 1 .AND. z ind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )140 DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 141 141 i_jpj = i_jpj - 1 142 142 END DO … … 221 221 ! 222 222 CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 223 CALL wrk_dealloc( jpj, z ind, zmsk )223 CALL wrk_dealloc( jpj, zswitch, zmsk ) 224 224 ! 225 225 IF( nn_timing == 1 ) CALL timing_stop('limdyn') … … 241 241 !!------------------------------------------------------------------- 242 242 INTEGER :: ios ! Local integer output status for namelist read 243 NAMELIST/namicedyn/ epsd, om, cw, angvg,pstar, &243 NAMELIST/namicedyn/ epsd, om, cw, pstar, & 244 244 & c_rhg, creepl, ecc, ahi0, & 245 245 & nevp, relast, alphaevp, hminrhg … … 262 262 WRITE(numout,*) ' relaxation constant om = ', om 263 263 WRITE(numout,*) ' drag coefficient for oceanic stress cw = ', cw 264 WRITE(numout,*) ' turning angle for oceanic stress angvg = ', angvg265 264 WRITE(numout,*) ' first bulk-rheology parameter pstar = ', pstar 266 265 WRITE(numout,*) ' second bulk-rhelogy parameter c_rhg = ', c_rhg … … 274 273 ENDIF 275 274 ! 276 IF( angvg /= 0._wp ) THEN277 CALL ctl_warn( 'lim_dyn_init: turning angle for oceanic stress not properly coded for EVP ', &278 & '(see limsbc module). We force angvg = 0._wp' )279 angvg = 0._wp280 ENDIF281 282 275 usecc2 = 1._wp / ( ecc * ecc ) 283 276 rhoco = rau0 * cw 284 angvg = angvg * rad285 sangvg = SIN( angvg )286 cangvg = COS( angvg )287 pstarh = pstar * 0.5_wp288 277 289 278 ! elastic damping -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r4924 r4972 43 43 PUBLIC lim_itd_me_alloc ! called by iceini.F90 44 44 45 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values46 REAL(wp) :: epsi10 = 1.e-10_wp ! constant values47 REAL(wp) :: epsi06 = 1.e-06_wp ! constant values48 49 45 !----------------------------------------------------------------------- 50 46 ! Variables shared among ridging subroutines … … 851 847 INTEGER :: ij ! horizontal index, combines i and j loops 852 848 INTEGER :: icells ! number of cells with aicen > puny 853 REAL(wp) :: zindb ! local scalar854 849 REAL(wp) :: hL, hR, farea, zdummy, zdummy0, ztmelts ! left and right limits of integration 855 850 REAL(wp) :: zsstK ! SST in Kelvin -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4924 r4972 46 46 PUBLIC lim_itd_shiftice 47 47 48 REAL(wp) :: epsi10 = 1.e-10_wp !49 REAL(wp) :: epsi06 = 1.e-6_wp !50 51 48 !!---------------------------------------------------------------------- 52 49 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) … … 156 153 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 157 154 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 158 REAL(wp) :: zx3, zareamin , zindb! - -155 REAL(wp) :: zx3, zareamin ! - - 159 156 CHARACTER (len = 15) :: fieldid 160 157 … … 219 216 DO jj = 1, jpj 220 217 DO ji = 1, jpi 221 zindb= 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes222 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb223 zindb= 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes224 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * zindb218 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 219 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 220 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 221 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 225 222 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 226 223 END DO … … 589 586 REAL(wp) :: zdo_aice ! ice age times volume transferred 590 587 REAL(wp) :: zdaTsf ! aicen*Tsfcn transferred 591 REAL(wp) :: zindsn ! snow or not592 REAL(wp) :: zindb ! ice or not593 588 594 589 INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions … … 717 712 718 713 jl1 = zdonor(ii,ij,jl) 719 zindb= MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) )720 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * zindb714 rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 715 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * rswitch 721 716 IF( jl1 == jl) THEN ; jl2 = jl1+1 722 717 ELSE ; jl2 = jl … … 814 809 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 815 810 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 816 zindsn= 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes811 rswitch = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 817 812 ELSE 818 813 ht_i(ji,jj,jl) = 0._wp -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r4924 r4972 50 50 PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2) 51 51 52 REAL(wp) :: epsi10 = 1.e-10_wp !53 54 52 !! * Substitutions 55 53 # include "vectopt_loop_substitute.h90" … … 119 117 CHARACTER (len=50) :: charout 120 118 REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta ! 121 REAL(wp) :: za, zstms, zsang, zmask ! local scalars 119 REAL(wp) :: za, zstms, zmask ! local scalars 120 REAL(wp) :: zc1, zc2, zc3 ! ice mass 122 121 123 122 REAL(wp) :: dtevp ! time step for subcycling … … 125 124 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars 126 125 REAL(wp) :: zu_ice2, zv_ice1 ! 127 REAL(wp) :: zddc, zdtc, zzdst ! delta on corners and on centre 126 REAL(wp) :: zddc, zdtc ! delta on corners and on centre 127 REAL(wp) :: zdst ! shear at the center of the grid point 128 128 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface 129 129 REAL(wp) :: sigma1, sigma2 ! internal ice stress 130 130 131 131 REAL(wp) :: zresm ! Maximal error on ice velocity 132 REAL(wp) :: zindb ! ice (1) or not (0)133 132 REAL(wp) :: zdummy ! dummy argument 134 133 REAL(wp) :: zintb, zintn ! dummy argument … … 140 139 REAL(wp), POINTER, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points 141 140 REAL(wp), POINTER, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays 142 REAL(wp), POINTER, DIMENSION(:,:) :: zc1 ! ice mass143 REAL(wp), POINTER, DIMENSION(:,:) :: zusw ! temporary weight for ice strength computation144 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce1, v_oce1 ! ocean u/v component on U points 145 142 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2, v_oce2 ! ocean u/v component on V points … … 147 144 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses 148 145 149 REAL(wp), POINTER, DIMENSION(:,:) :: zd d , zdt ! Divergence andtension at centre of grid cells146 REAL(wp), POINTER, DIMENSION(:,:) :: zdt ! tension at centre of grid cells 150 147 REAL(wp), POINTER, DIMENSION(:,:) :: zds ! Shear on northeast corner of grid cells 151 REAL(wp), POINTER, DIMENSION(:,:) :: zdst ! Shear on centre of grid cells152 REAL(wp), POINTER, DIMENSION(:,:) :: deltat, deltac ! Delta at centre and corners of grid cells153 148 REAL(wp), POINTER, DIMENSION(:,:) :: zs1 , zs2 ! Diagonal stress tensor components zs1 and zs2 154 149 REAL(wp), POINTER, DIMENSION(:,:) :: zs12 ! Non-diagonal stress tensor component zs12 … … 160 155 161 156 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 162 CALL wrk_alloc( jpi,jpj, zc1 , u_oce1, u_oce2, u_ice2, zusw, v_oce1 , v_oce2, v_ice1 )163 CALL wrk_alloc( jpi,jpj, zf1 , deltat, zu_ice, zf2 , deltac, zv_ice , zdd , zdt , zds , zdst)164 CALL wrk_alloc( jpi,jpj, zd d , zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )157 CALL wrk_alloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1 ) 158 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 159 CALL wrk_alloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) 165 160 166 161 #if defined key_lim2 && ! defined key_lim2_vp … … 179 174 ! 180 175 !------------------------------------------------------------------------------! 181 ! 1) Ice -Snow mass (zc1), icestrength (zpresh) !176 ! 1) Ice strength (zpresh) ! 182 177 !------------------------------------------------------------------------------! 183 178 ! 184 179 ! Put every vector to 0 185 zpresh (:,:) = 0._wp ; zc1 (:,:) = 0._wp 180 delta_i(:,:) = 0._wp ; 181 zpresh (:,:) = 0._wp ; 186 182 zpreshc(:,:) = 0._wp 187 183 u_ice2 (:,:) = 0._wp ; v_ice1(:,:) = 0._wp 188 zdd (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp 184 divu_i (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp 185 shear_i(:,:) = 0._wp 189 186 190 187 #if defined key_lim3 … … 196 193 !CDIR NOVERRCHK 197 194 DO ji = 1 , jpi 198 zc1(ji,jj) = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )199 195 #if defined key_lim3 200 196 zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj) … … 218 214 & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 219 215 & tms(ji,jj) * wght(ji+1,jj+1,1,1) 220 zusw(ji,jj) = 1.0 / MAX( zstms, epsd )221 216 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 222 217 & zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 223 218 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 224 219 & zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 225 & ) * zusw(ji,jj)220 & ) / MAX( zstms, epsd ) 226 221 END DO 227 222 END DO … … 265 260 DO ji = fs_2, fs_jpim1 266 261 262 zc1 = tms(ji ,jj ) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 263 zc2 = tms(ji+1,jj ) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 264 zc3 = tms(ji ,jj+1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 265 267 266 zt11 = tms(ji ,jj) * e1t(ji ,jj) 268 267 zt12 = tms(ji+1,jj) * e1t(ji+1,jj) … … 275 274 276 275 ! Mass, coriolis coeff. and currents 277 zmass1(ji,jj) = ( zt12*zc1 (ji,jj) + zt11*zc1(ji+1,jj)) / (zt11+zt12+epsd)278 zmass2(ji,jj) = ( zt22*zc1 (ji,jj) + zt21*zc1(ji,jj+1)) / (zt21+zt22+epsd)276 zmass1(ji,jj) = ( zt12*zc1 + zt11*zc2 ) / (zt11+zt12+epsd) 277 zmass2(ji,jj) = ( zt22*zc1 + zt21*zc3 ) / (zt21+zt22+epsd) 279 278 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) ) & 280 279 & / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd ) … … 344 343 ! 345 344 !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 346 !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells345 !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells 347 346 !- zds(:,:): shear on northeast corner of grid cells 348 347 ! … … 353 352 ! bugs (Martin, for Miguel). 354 353 ! 355 !- ALSO: arrays zd d, zdt, zds and delta could354 !- ALSO: arrays zdt, zds and delta could 356 355 ! be removed in the future to minimise memory demand. 357 356 ! … … 361 360 ! 362 361 ! 363 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) &364 & -e2u(ji-1,jj)*u_ice(ji-1,jj) &365 & +e1v(ji,jj)*v_ice(ji,jj) &366 & -e1v(ji,jj-1)*v_ice(ji,jj-1) &367 & ) &368 & / area(ji,jj)362 divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 363 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 364 & +e1v(ji,jj)*v_ice(ji,jj) & 365 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 366 & ) & 367 & / area(ji,jj) 369 368 370 369 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & … … 408 407 409 408 !- Calculate Delta at centre of grid cells 410 z zdst = ( e2u(ji , jj) * v_ice1(ji ,jj) &409 zdst = ( e2u(ji , jj) * v_ice1(ji ,jj) & 411 410 & - e2u(ji-1, jj) * v_ice1(ji-1,jj) & 412 411 & + e1v(ji, jj ) * u_ice2(ji,jj ) & … … 415 414 & / area(ji,jj) 416 415 417 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 ) 418 ! MV rewriting 419 ! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 420 !!gm faster to replace the line above with simply: 421 !! deltat(ji,jj) = MAX( delta, creepl ) 422 !!gm end 423 deltat(ji,jj) = delta + creepl 424 ! END MV 416 delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 417 delta_i(ji,jj) = delta + creepl 425 418 !-Calculate stress tensor components zs1 and zs2 426 419 !-at centre of grid cells (see section 3.5 of CICE user's guide). 427 !zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) + & 428 ! & ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) ) & 429 ! & / ( 1._wp + alphaevp * dtotel ) 430 431 !zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) - & 432 ! zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) ) & 433 ! & / ( 1._wp + alphaevp * ecc2 * dtotel ) 434 435 ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 436 zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) ) & 420 zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) ) & 437 421 & * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 438 zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta t(ji,jj) * zpresh(ji,jj) ) ) &422 zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) ) & 439 423 & / ( 1._wp + dtotel ) 440 424 … … 468 452 & / ( e1f(ji,jj) * e2f(ji,jj) ) 469 453 470 deltac(ji,jj)= SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl454 zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 471 455 472 456 !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 473 !zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) / &474 ! & ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) ) &475 ! & / ( 1._wp + alphaevp * ecc2 * dtotel )476 477 ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)478 457 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * & 479 & ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj)) * zpreshc(ji,jj) ) ) &458 & ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) ) & 480 459 & / ( 1.0 + dtotel ) 481 460 … … 513 492 DO ji = fs_2, fs_jpim1 514 493 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 515 zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg516 494 z0 = zmass1(ji,jj)/dtevp 517 495 … … 523 501 (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 524 502 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 525 za*( cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))526 zcca = z0+za *cangvg527 zccb = zcorl1(ji,jj) +za*zsang503 za*(u_oce1(ji,jj)) 504 zcca = z0+za 505 zccb = zcorl1(ji,jj) 528 506 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 529 507 … … 536 514 #endif 537 515 #if defined key_bdy 538 ! clem: change u_ice and v_ice at the boundary for each iteration539 516 CALL bdy_ice_lim_dyn( 'U' ) 540 517 #endif … … 546 523 547 524 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 548 zsang = SIGN(1.0,fcor(ji,jj))*sangvg549 525 z0 = zmass2(ji,jj)/dtevp 550 526 ! SB modif because ocean has no slip boundary condition … … 555 531 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 556 532 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 557 za2ct(ji,jj) + za*( cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))558 zcca = z0+za *cangvg559 zccb = zcorl2(ji,jj) +za*zsang533 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 534 zcca = z0+za 535 zccb = zcorl2(ji,jj) 560 536 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 561 537 … … 568 544 #endif 569 545 #if defined key_bdy 570 ! clem: change u_ice and v_ice at the boundary for each iteration571 546 CALL bdy_ice_lim_dyn( 'V' ) 572 547 #endif … … 578 553 DO ji = fs_2, fs_jpim1 579 554 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 580 zsang = SIGN(1.0,fcor(ji,jj))*sangvg581 555 z0 = zmass2(ji,jj)/dtevp 582 556 ! SB modif because ocean has no slip boundary condition … … 588 562 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 589 563 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 590 za2ct(ji,jj) + za*( cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))591 zcca = z0+za *cangvg592 zccb = zcorl2(ji,jj) +za*zsang564 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 565 zcca = z0+za 566 zccb = zcorl2(ji,jj) 593 567 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 594 568 … … 601 575 #endif 602 576 #if defined key_bdy 603 ! clem: change u_ice and v_ice at the boundary for each iteration604 577 CALL bdy_ice_lim_dyn( 'V' ) 605 578 #endif … … 610 583 DO ji = fs_2, fs_jpim1 611 584 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 612 zsang = SIGN(1.0,fcor(ji,jj))*sangvg613 585 z0 = zmass1(ji,jj)/dtevp 614 ! SB modif because ocean has no slip boundary condition615 ! GG Bug616 ! zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) &617 ! & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &618 ! & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)619 586 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 620 587 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & … … 624 591 (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 625 592 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 626 za*( cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))627 zcca = z0+za *cangvg628 zccb = zcorl1(ji,jj) +za*zsang593 za*(u_oce1(ji,jj)) 594 zcca = z0+za 595 zccb = zcorl1(ji,jj) 629 596 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 630 597 END DO ! ji … … 636 603 #endif 637 604 #if defined key_bdy 638 ! clem: change u_ice and v_ice at the boundary for each iteration639 605 CALL bdy_ice_lim_dyn( 'U' ) 640 606 #endif … … 666 632 !CDIR NOVERRCHK 667 633 DO ji = fs_2, fs_jpim1 668 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )669 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )670 634 zdummy = vt_i(ji,jj) 671 635 IF ( zdummy .LE. hminrhg ) THEN … … 683 647 #endif 684 648 #if defined key_bdy 685 ! clem: change u_ice and v_ice at the boundary686 649 CALL bdy_ice_lim_dyn( 'U' ) 687 650 CALL bdy_ice_lim_dyn( 'V' ) … … 690 653 DO jj = k_j1+1, k_jpj-1 691 654 DO ji = fs_2, fs_jpim1 692 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )693 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )694 655 zdummy = vt_i(ji,jj) 695 656 IF ( zdummy .LE. hminrhg ) THEN … … 713 674 !CDIR NOVERRCHK 714 675 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi 715 !- zdd(:,:), zdt(:,:): divergence and tension at centre676 !- divu_i(:,:), zdt(:,:): divergence and tension at centre 716 677 !- zds(:,:): shear on northeast corner of grid cells 717 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )718 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )719 678 zdummy = vt_i(ji,jj) 720 679 IF ( zdummy .LE. hminrhg ) THEN 721 680 722 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) &723 & -e2u(ji-1,jj)*u_ice(ji-1,jj) &724 & +e1v(ji,jj)*v_ice(ji,jj) &725 & -e1v(ji,jj-1)*v_ice(ji,jj-1) &726 & ) &727 & / area(ji,jj)681 divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 682 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 683 & +e1v(ji,jj)*v_ice(ji,jj) & 684 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 685 & ) & 686 & / area(ji,jj) 728 687 729 688 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & … … 747 706 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 748 707 749 zdst (ji,jj)= ( e2u( ji , jj ) * v_ice1(ji ,jj ) &708 zdst = ( e2u( ji , jj ) * v_ice1(ji ,jj ) & 750 709 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj ) & 751 710 & + e1v( ji , jj ) * u_ice2(ji ,jj ) & 752 711 & - e1v( ji , jj-1 ) * u_ice2(ji ,jj-1) ) / area(ji,jj) 753 712 754 ! deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) & 755 ! & + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 756 ! & ) + creepl 757 ! MV rewriting 758 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 ) 759 deltat(ji,jj) = delta + creepl 760 ! END MV 713 delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 714 delta_i(ji,jj) = delta + creepl 761 715 762 716 ENDIF ! zdummy … … 773 727 DO jj = k_j1+1, k_jpj-1 774 728 DO ji = fs_2, fs_jpim1 775 divu_i (ji,jj) = zdd (ji,jj)776 delta_i(ji,jj) = deltat(ji,jj)777 729 ! begin TECLIM change 778 zdst (ji,jj)= ( e2u( ji , jj ) * v_ice1(ji,jj) &730 zdst= ( e2u( ji , jj ) * v_ice1(ji,jj) & 779 731 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 780 732 & + e1v( ji , jj ) * u_ice2(ji,jj) & 781 733 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 782 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst (ji,jj) * zdst(ji,jj))734 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 783 735 ! end TECLIM change 784 736 END DO … … 834 786 ! 835 787 CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 836 CALL wrk_dealloc( jpi,jpj, zc1 , u_oce1, u_oce2, u_ice2, zusw, v_oce1 , v_oce2, v_ice1 )837 CALL wrk_dealloc( jpi,jpj, zf1 , deltat, zu_ice, zf2 , deltac, zv_ice , zdd , zdt , zds , zdst)838 CALL wrk_dealloc( jpi,jpj, zd d , zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )788 CALL wrk_dealloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1 ) 789 CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 790 CALL wrk_dealloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) 839 791 840 792 END SUBROUTINE lim_rhg -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
r4924 r4972 308 308 INTEGER :: ji, jj, jk, jl, indx 309 309 REAL(wp) :: zfice, ziter 310 REAL(wp) :: zs_inf, z_slope_s, zsmax, zsmin, zalpha , zindb! local scalars used for the salinity profile311 REAL(wp), POINTER, DIMENSION(:) :: zs_zero310 REAL(wp) :: zs_inf, z_slope_s, zsmax, zsmin, zalpha ! local scalars used for the salinity profile 311 REAL(wp), POINTER, DIMENSION(:) :: zs_zero 312 312 REAL(wp), POINTER, DIMENSION(:,:) :: z2d 313 313 CHARACTER(len=15) :: znam -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r4967 r4972 50 50 PUBLIC lim_sbc_flx ! called by sbc_ice_lim 51 51 PUBLIC lim_sbc_tau ! called by sbc_ice_lim 52 53 REAL(wp) :: epsi10 = 1.e-10 !54 REAL(wp) :: epsi20 = 1.e-20 !55 52 56 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_oce, vtau_oce ! air-ocean surface i- & j-stress [N/m2] … … 108 105 INTEGER :: ji, jj, jl, jk ! dummy loop indices 109 106 ! 110 REAL(wp) :: z inda, zemp! local scalars107 REAL(wp) :: zemp ! local scalars 111 108 REAL(wp) :: zf_mass ! Heat flux associated with mass exchange ice->ocean (W.m-2) 112 109 REAL(wp) :: zfcm1 ! New solar flux received by the ocean … … 131 128 ! heat flux at the ocean surface ! 132 129 !------------------------------------------! 133 zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( 1._wp - pfrld(ji,jj) ) ) ) ! 1 if ice134 135 130 ! Solar heat flux reaching the ocean = zfcm1 (W.m-2) 136 131 !--------------------------------------------------- -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4946 r4972 50 50 PUBLIC lim_thd ! called by limstp module 51 51 PUBLIC lim_thd_init ! called by iceini module 52 53 REAL(wp) :: epsi10 = 1.e-10_wp !54 52 55 53 !! * Substitutions … … 89 87 REAL(wp) :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 90 88 REAL(wp) :: zch = 0.0057_wp ! heat transfer coefficient 91 REAL(wp) :: z inda, zindb, zareamin89 REAL(wp) :: zareamin 92 90 REAL(wp) :: zfric_u, zqld, zqfr 93 91 ! … … 103 101 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 104 102 105 !------------------------------------------------------------------------------! 106 ! 1) Initialization of diagnostic variables ! 107 !------------------------------------------------------------------------------! 103 !------------------------------------------------------------------------! 104 ! 1) Initialization of some variables ! 105 !------------------------------------------------------------------------! 106 ftr_ice(:,:,:) = 0._wp ! part of solar radiation transmitted through the ice 107 108 108 109 109 !-------------------- … … 116 116 DO ji = 1, jpi 117 117 !0 if no ice and 1 if yes 118 zindb= 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )118 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) 119 119 !Energy of melting q(S,T) [J.m-3] 120 e_i(ji,jj,jk,jl) = zindb* e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i )120 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 121 121 !convert units ! very important that this line is here 122 122 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac … … 128 128 DO ji = 1, jpi 129 129 !0 if no ice and 1 if yes 130 zindb= 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) )130 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) ) 131 131 !Energy of melting q(S,T) [J.m-3] 132 e_s(ji,jj,jk,jl) = zindb* e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s )132 e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 133 133 !convert units ! very important that this line is here 134 134 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac … … 165 165 !CDIR NOVERRCHK 166 166 DO ji = 1, jpi 167 zinda= tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice167 rswitch = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 168 168 ! 169 169 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 201 201 fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90 202 202 qlead(ji,jj) = 0._wp 203 ELSE 204 fhld (ji,jj) = 0._wp 203 205 ENDIF 204 206 ! … … 206 208 !clem zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 207 209 zfric_u = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 208 fhtur(ji,jj) = MAX( 0._wp, zinda* rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2210 fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2 209 211 ! upper bound for fhtur: we do not want SST to drop below Tfreeze. 210 212 ! So we say that the heat retrieved from the ocean (fhtur+fhld) must be < to the heat necessary to reach Tfreeze (zqfr) 211 213 ! This is not a clean budget, so that should be corrected at some point 212 fhtur(ji,jj) = zinda* MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )214 fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 213 215 214 216 ! ----------------------------------------- … … 418 420 CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj ) 419 421 CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj ) 420 !421 IF( num_sal == 2 ) THEN422 422 CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj ) 423 ENDIF424 423 425 424 CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj ) … … 436 435 CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb) , jpi, jpj ) 437 436 ! 438 !+++++ temporary stuff for a dummy version439 CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj )440 CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj )441 CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new (1:nbpb) , jpi, jpj )442 CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0 (1:nbpb) , jpi, jpj )443 !+++++444 437 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 445 438 CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) … … 536 529 !! 537 530 INTEGER :: ji, jk ! dummy loop indices 538 REAL(wp) :: ztmelts, z switch, zaaa, zbbb, zccc, zdiscrim ! local scalar531 REAL(wp) :: ztmelts, zaaa, zbbb, zccc, zdiscrim ! local scalar 539 532 !!------------------------------------------------------------------- 540 533 ! Recover ice temperature … … 547 540 zccc = lfus * ( ztmelts - rtt ) 548 541 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 549 t_i_1d(ji,jk) 542 t_i_1d(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 550 543 551 544 ! mask temperature 552 zswitch= 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )553 t_i_1d(ji,jk) = zswitch * t_i_1d(ji,jk) + ( 1._wp - zswitch ) * rtt545 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 546 t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rtt 554 547 END DO 555 548 END DO -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4946 r4972 32 32 PUBLIC lim_thd_dh ! called by lim_thd 33 33 34 REAL(wp) :: epsi20 = 1.e-20 ! constant values35 REAL(wp) :: epsi10 = 1.e-10 !36 37 34 !!---------------------------------------------------------------------- 38 35 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) … … 111 108 112 109 ! mass and salt flux (clem) 113 REAL(wp) :: zdvres, zswitch_sal , zswitch110 REAL(wp) :: zdvres, zswitch_sal 114 111 115 112 ! Heat conservation 116 113 INTEGER :: num_iter_max 117 REAL(wp) :: zinda, zindq, zindh118 REAL(wp), POINTER, DIMENSION(:) :: zintermelt ! debug119 114 120 115 !!------------------------------------------------------------------ … … 128 123 CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 129 124 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 130 CALL wrk_alloc( jpij, zintermelt )131 125 CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) 132 126 CALL wrk_alloc( jpij, icount ) … … 147 141 zh_i (:,:) = 0._wp 148 142 zdeltah (:,:) = 0._wp 149 zintermelt(:) = 0._wp150 143 icount (:) = 0 151 144 … … 165 158 ! 166 159 DO ji = kideb, kiut 167 zinda= 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )168 ztmelts = zinda * rtt + ( 1._wp - zinda) * rtt160 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 161 ztmelts = rswitch * rtt + ( 1._wp - rswitch ) * rtt 169 162 170 163 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) … … 254 247 ! thickness change 255 248 IF( zdh_s_pre(ji) > 0._wp ) THEN 256 zindq= 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) )257 zdh_s_mel (ji) = - zindq* zq_su(ji) / MAX( zqprec(ji) , epsi20 )249 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) ) 250 zdh_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 258 251 zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting 259 252 ! heat used to melt snow (W.m-2, >0) … … 275 268 DO ji = kideb, kiut 276 269 ! thickness change 277 zindh= 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )278 zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20) )279 zdeltah (ji,jk) = - zindh * zindq* zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )270 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 271 rswitch = rswitch * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) ) ) 272 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 280 273 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting 281 274 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) … … 331 324 DO jk = 1, nlay_s 332 325 DO ji = kideb,kiut 333 zindh= MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 ) )334 q_s_1d(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_1d(ji), epsi20 ) * &326 rswitch = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 ) ) 327 q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) * & 335 328 & ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + & 336 329 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) ) … … 382 375 ! => icount=0 : no layer has vanished 383 376 ! => icount=5 : 5 layers have vanished 384 zindh = NINT( MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk)) ) ) )385 icount(ji) = icount(ji) + zindh377 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 378 icount(ji) = icount(ji) + NINT( rswitch ) 386 379 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 387 380 … … 515 508 516 509 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 517 zintermelt(ji) = 1._wp518 510 519 511 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) … … 602 594 ! 603 595 ! ! excessive energy is sent to lateral ablation 604 ! zinda= MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) )605 ! zq_1cat(ji) = zinda* rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0596 ! rswitch = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) ) 597 ! zq_1cat(ji) = rswitch * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0 606 598 ! 607 599 ! ! correct salt and mass fluxes … … 668 660 ! new salinity difference stored (to be used in limthd_ent.F90) 669 661 IF ( num_sal == 2 ) THEN 670 zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) )662 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 671 663 ! salinity dif due to snow-ice formation 672 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch664 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch 673 665 ! salinity dif due to bottom growth 674 666 IF ( zf_tt(ji) < 0._wp ) THEN 675 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch667 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch 676 668 ENDIF 677 669 ENDIF … … 710 702 !clem bug: we should take snow into account here 711 703 DO ji = kideb, kiut 712 zindh= 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )713 t_su_1d(ji) = zindh * t_su_1d(ji) + ( 1.0 - zindh ) * rtt704 rswitch = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 705 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rtt 714 706 END DO ! ji 715 707 … … 717 709 DO ji = kideb,kiut 718 710 ! mask enthalpy 719 zinda= MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) )720 q_s_1d(ji,jk) = ( 1.0 - zinda) * q_s_1d(ji,jk)711 rswitch = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) 712 q_s_1d(ji,jk) = ( 1.0 - rswitch ) * q_s_1d(ji,jk) 721 713 ! recalculate t_s_1d from q_s_1d 722 t_s_1d(ji,jk) = rtt + ( 1._wp - zinda) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )714 t_s_1d(ji,jk) = rtt + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 723 715 END DO 724 716 END DO … … 726 718 CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 727 719 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 728 CALL wrk_dealloc( jpij, zintermelt )729 720 CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) 730 721 CALL wrk_dealloc( jpij, icount ) -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4964 r4972 32 32 PUBLIC lim_thd_dif ! called by lim_thd 33 33 34 REAL(wp) :: epsi10 = 1.e-10_wp !35 34 !!---------------------------------------------------------------------- 36 35 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 141 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 142 141 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: z indterm ! Independent term144 REAL(wp), POINTER, DIMENSION(:,:) :: z indtbis ! temporary independent term142 REAL(wp), POINTER, DIMENSION(:,:) :: zswiterm ! Independent term 143 REAL(wp), POINTER, DIMENSION(:,:) :: zswitbis ! temporary independent term 145 144 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis 146 145 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms … … 154 153 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 155 154 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 156 CALL wrk_alloc( jpij, nlay_i+3, z indterm, zindtbis, zdiagbis )155 CALL wrk_alloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis ) 157 156 CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 158 157 … … 438 437 ztrid(ji,numeq,2) = 0. 439 438 ztrid(ji,numeq,3) = 0. 440 z indterm(ji,numeq)= 0.441 z indtbis(ji,numeq)= 0.439 zswiterm(ji,numeq)= 0. 440 zswitbis(ji,numeq)= 0. 442 441 zdiagbis(ji,numeq)= 0. 443 442 ENDDO … … 451 450 zkappa_i(ji,jk)) 452 451 ztrid(ji,numeq,3) = - zeta_i(ji,jk)*zkappa_i(ji,jk) 453 z indterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* &452 zswiterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* & 454 453 zradab_i(ji,jk) 455 454 END DO … … 463 462 + zkappa_i(ji,nlay_i-1) ) 464 463 ztrid(ji,numeq,3) = 0.0 465 z indterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* &464 zswiterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 466 465 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 467 466 * t_bo_1d(ji) ) … … 483 482 zkappa_s(ji,jk) ) 484 483 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 485 z indterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* &484 zswiterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* & 486 485 zradab_s(ji,jk) 487 486 END DO … … 490 489 IF ( nlay_i.eq.1 ) THEN 491 490 ztrid(ji,nlay_s+2,3) = 0.0 492 z indterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &491 zswiterm(ji,nlay_s+2) = zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 493 492 t_bo_1d(ji) 494 493 ENDIF … … 507 506 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 508 507 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 509 z indterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji)508 zswiterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji) 510 509 511 510 !!first layer of snow equation … … 513 512 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 514 513 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 515 z indterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)514 zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 516 515 517 516 ELSE … … 530 529 zkappa_s(ji,0) * zg1s ) 531 530 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 532 z indterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * &531 zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 533 532 ( zradab_s(ji,1) + & 534 533 zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) … … 554 553 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 555 554 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 556 z indterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji)555 zswiterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 557 556 558 557 !!first layer of ice equation … … 561 560 + zkappa_i(ji,0) * zg1 ) 562 561 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 563 z indterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)562 zswiterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 564 563 565 564 !!case of only one layer in the ice (surface & ice equations are altered) … … 574 573 ztrid(ji,numeqmin(ji)+1,3) = 0.0 575 574 576 z indterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* &575 zswiterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* & 577 576 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 578 577 ENDIF … … 594 593 zg1) 595 594 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 596 z indterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &595 zswiterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 597 596 zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 598 597 … … 603 602 zkappa_i(ji,1)) 604 603 ztrid(ji,numeqmin(ji),3) = 0.0 605 z indterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* &604 zswiterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* & 606 605 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 607 606 + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 … … 627 626 628 627 DO ji = kideb , kiut 629 z indtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji))628 zswitbis(ji,numeqmin(ji)) = zswiterm(ji,numeqmin(ji)) 630 629 zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2) 631 630 minnumeqmin = MIN(numeqmin(ji),minnumeqmin) … … 638 637 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 639 638 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 640 z indtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* &641 z indtbis(ji,numeq-1)/zdiagbis(ji,numeq-1)639 zswitbis(ji,numeq) = zswiterm(ji,numeq) - ztrid(ji,numeq,1)* & 640 zswitbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 642 641 END DO 643 642 END DO … … 645 644 DO ji = kideb , kiut 646 645 ! ice temperatures 647 t_i_1d(ji,nlay_i) = z indtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))646 t_i_1d(ji,nlay_i) = zswitbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 648 647 END DO 649 648 … … 651 650 DO ji = kideb , kiut 652 651 jk = numeq - nlay_s - 1 653 t_i_1d(ji,jk) = (z indtbis(ji,numeq) - ztrid(ji,numeq,3)* &652 t_i_1d(ji,jk) = (zswitbis(ji,numeq) - ztrid(ji,numeq,3)* & 654 653 t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 655 654 END DO … … 659 658 ! snow temperatures 660 659 IF (ht_s_1d(ji).GT.0._wp) & 661 t_s_1d(ji,nlay_s) = (z indtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &660 t_s_1d(ji,nlay_s) = (zswitbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 662 661 * t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 663 662 * MAX(0.0,SIGN(1.0,ht_s_1d(ji))) … … 667 666 ztsubit(ji) = t_su_1d(ji) 668 667 IF( t_su_1d(ji) < ztfs(ji) ) & 669 t_su_1d(ji) = ( z indtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) &668 t_su_1d(ji) = ( zswitbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) & 670 669 & + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 671 670 END DO … … 780 779 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 781 780 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 782 CALL wrk_dealloc( jpij, nlay_i+3, z indterm, zindtbis, zdiagbis )781 CALL wrk_dealloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis ) 783 782 CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 784 783 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) … … 797 796 ! 798 797 INTEGER :: ji, jk ! dummy loop indices 799 REAL(wp) :: ztmelts , zindb! local scalar798 REAL(wp) :: ztmelts ! local scalar 800 799 !!------------------------------------------------------------------- 801 800 ! … … 803 802 DO ji = kideb, kiut 804 803 ztmelts = - tmut * s_i_1d(ji,jk) + rtt 805 zindb= MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) )804 rswitch = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) ) 806 805 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 807 & + lfus * ( 1.0 - zindb* ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) ) &806 & + lfus * ( 1.0 - rswitch * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) ) & 808 807 & - rcp * ( ztmelts-rtt ) ) 809 808 END DO -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r4924 r4972 38 38 PUBLIC lim_thd_ent ! called by limthd and limthd_lac 39 39 40 REAL(wp) :: epsi20 = 1.e-20 ! constant values41 REAL(wp) :: epsi10 = 1.e-10 ! constant values42 43 40 !!---------------------------------------------------------------------- 44 41 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 79 76 INTEGER :: ji ! dummy loop indices 80 77 INTEGER :: jk0, jk1 ! old/new layer indices 81 REAL(wp) :: zswitch82 78 ! 83 79 REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces … … 137 133 DO jk1 = 1, nlay_i 138 134 DO ji = kideb, kiut 139 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) )140 qnew(ji,jk1) = zswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi10 )135 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) ) 136 qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi10 ) 141 137 ENDDO 142 138 ENDDO -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4946 r4972 37 37 38 38 PUBLIC lim_thd_lac ! called by lim_thd 39 40 REAL(wp) :: epsi10 = 1.e-10_wp !41 REAL(wp) :: epsi20 = 1.e-20_wp !42 39 43 40 !!---------------------------------------------------------------------- … … 77 74 INTEGER :: nbpac ! local integers 78 75 INTEGER :: ii, ij, iter ! - - 79 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, z indb, zinda, zde! local scalars76 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zde ! local scalars 80 77 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 81 78 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - … … 132 129 DO ji = 1, jpi 133 130 !Energy of melting q(S,T) [J.m-3] 134 zindb= 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes135 e_i(ji,jj,jk,jl) = zindb* e_i(ji,jj,jk,jl) &131 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 132 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) & 136 133 & / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i, wp ) 137 134 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac … … 189 186 ! Frazil ice velocity 190 187 !--------------------- 191 zindb= MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )192 zvfrx = zindb* zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )193 zvfry = zindb* zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )188 rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 189 zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 190 zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 194 191 195 192 !------------------- … … 197 194 !------------------- 198 195 ! C-grid ice velocity 199 zindb = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 200 zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) & 201 & + u_ice(ji,jj ) * tmu(ji ,jj ) ) * 0.5_wp 202 zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) & 203 & + v_ice(ji,jj ) * tmv(ji ,jj ) ) * 0.5_wp 196 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 197 zvgx = rswitch * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) + u_ice(ji,jj) * tmu(ji,jj) ) * 0.5_wp 198 zvgy = rswitch * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) + v_ice(ji,jj) * tmv(ji,jj) ) * 0.5_wp 204 199 205 200 !----------------------------------- … … 391 386 392 387 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 393 zinda= 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )394 zfrazb = zinda* ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb388 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 389 zfrazb = rswitch * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 395 390 zv_frazb(ji) = zfrazb * zv_newice(ji) 396 391 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) … … 447 442 DO ji = 1, nbpac 448 443 jl = jcat(ji) 449 zinda= MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )444 rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 450 445 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 451 446 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) ) & 452 & * zinda/ MAX( zv_i_1d(ji,jl), epsi20 )447 & * rswitch / MAX( zv_i_1d(ji,jl), epsi20 ) 453 448 END DO 454 449 END DO … … 471 466 ! new volumes including lateral/bottom accretion + residual 472 467 DO ji = 1, nbpac 473 zinda= MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )474 zv_newfra = zinda* ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )475 za_i_1d(ji,jl) = zinda* za_i_1d(ji,jl)468 rswitch = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 469 zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 470 za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl) 476 471 zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 477 472 ! for remapping … … 488 483 DO jl = 1, jpl 489 484 DO ji = 1, nbpac 490 zindb= 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes491 zoa_i_1d(ji,jl) = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb485 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes 486 zoa_i_1d(ji,jl) = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch 492 487 END DO 493 488 END DO -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4963 r4972 37 37 PUBLIC lim_trp ! called by ice_step 38 38 39 REAL(wp) :: epsi10 = 1.e-10_wp40 REAL(wp) :: epsi20 = 1.e-20_wp41 42 39 !! * Substitution 43 40 # include "vectopt_loop_substitute.h90" … … 65 62 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 66 63 INTEGER :: initad ! number of sub-timestep for the advection 67 INTEGER :: ierr ! error status 68 REAL(wp) :: zindb , zindsn , zindic, zindh, zinda ! local scalar 69 REAL(wp) :: zcfl , zusnit ! - - 70 REAL(wp) :: zsal , zage ! - - 64 REAL(wp) :: zcfl , zusnit ! - - 71 65 ! 72 66 REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow 73 67 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 74 68 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e 75 ! mass and salt flux (clem) 76 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold ! old ice volume... 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 78 REAL(wp), POINTER, DIMENSION(:,:) :: zeiold, zesold ! old enthalpies 69 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold ! old ice volume... 70 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 71 REAL(wp), POINTER, DIMENSION(:,:) :: zeiold, zesold ! old enthalpies 79 72 REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei 80 73 ! … … 341 334 DO jj = 1, jpj 342 335 DO ji = 1, jpi 343 zindb= MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) )336 rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 344 337 345 338 zvi = zs0ice(ji,jj,jl) … … 349 342 ! 350 343 ! Remove very small areas 351 v_s(ji,jj,jl) = zindb* zs0sn (ji,jj,jl)352 v_i(ji,jj,jl) = zindb* zs0ice(ji,jj,jl)353 a_i(ji,jj,jl) = zindb* zs0a (ji,jj,jl)354 e_s(ji,jj,1,jl) = zindb* zs0c0 (ji,jj,jl)344 v_s(ji,jj,jl) = rswitch * zs0sn (ji,jj,jl) 345 v_i(ji,jj,jl) = rswitch * zs0ice(ji,jj,jl) 346 a_i(ji,jj,jl) = rswitch * zs0a (ji,jj,jl) 347 e_s(ji,jj,1,jl) = rswitch * zs0c0 (ji,jj,jl) 355 348 ! Ice salinity and age 356 349 IF( num_sal == 2 ) THEN 357 350 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) ) 358 351 ENDIF 359 oa_i(ji,jj,jl) = MAX( zindb* zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl)352 oa_i(ji,jj,jl) = MAX( rswitch * zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl) 360 353 361 354 ! Update fluxes … … 372 365 DO jj = 1, jpj 373 366 DO ji = 1, jpi 374 zindb= MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) )367 rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 375 368 zei = zs0e(ji,jj,jk,jl) 376 e_i(ji,jj,jk,jl) = zindb* MAX( 0._wp, zs0e(ji,jj,jk,jl) )369 e_i(ji,jj,jk,jl) = rswitch * MAX( 0._wp, zs0e(ji,jj,jk,jl) ) 377 370 ! Update fluxes 378 371 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 … … 397 390 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 398 391 399 zindh = 1._wp392 rswitch = 1._wp 400 393 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 401 394 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 402 395 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 403 zindh =MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )404 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )396 rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 397 a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 405 398 ELSE 406 399 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 407 zindh =MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )408 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )400 rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 401 a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 409 402 ENDIF 410 403 411 ! small correction due to * zindh for a_i412 v_i (ji,jj,jl) = zindh * v_i (ji,jj,jl)413 v_s (ji,jj,jl) = zindh * v_s (ji,jj,jl)414 smv_i(ji,jj,jl) = zindh * smv_i(ji,jj,jl)415 e_s(ji,jj,1,jl) = zindh * e_s(ji,jj,1,jl)416 e_i(ji,jj,1:nlay_i,jl) = zindh * e_i(ji,jj,1:nlay_i,jl)404 ! small correction due to *rswitch for a_i 405 v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl) 406 v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl) 407 smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl) 408 e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl) 409 e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 417 410 418 411 ! Update mass fluxes … … 422 415 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 423 416 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 424 425 417 ENDIF 426 427 diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice428 diag_trp_vs(ji,jj) = diag_trp_vs(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * r1_rdtice429 430 418 END DO 431 419 END DO … … 438 426 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 439 427 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 440 END DO 441 END DO 442 443 ! --- agglomerate variables (clem) ----------------- 428 429 diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 430 diag_trp_vs(ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 431 END DO 432 END DO 433 434 ! --- agglomerate variables ----------------- 444 435 vt_i (:,:) = 0._wp 445 436 vt_s (:,:) = 0._wp … … 462 453 DO ji = 1, jpi 463 454 ! open water = 1 if at_i=0 464 zindb= MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )465 ato_i(ji,jj) = zindb + (1._wp - zindb) * zs0ow(ji,jj)455 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 456 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * zs0ow(ji,jj) 466 457 END DO 467 458 END DO -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r4924 r4972 50 50 PUBLIC lim_update1 ! routine called by ice_step 51 51 52 REAL(wp) :: epsi10 = 1.e-10_wp ! - -53 54 52 !! * Substitutions 55 53 # include "vectopt_loop_substitute.h90" -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r4924 r4972 47 47 PUBLIC lim_update2 ! routine called by ice_step 48 48 49 REAL(wp) :: epsi10 = 1.e-10_wp ! - -50 REAL(wp) :: epsi20 = 1.e-20_wp51 52 49 !! * Substitutions 53 50 # include "vectopt_loop_substitute.h90" -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r4924 r4972 66 66 PUBLIC lim_var_salprof1d ! 67 67 68 REAL(wp) :: epsi10 = 1.e-10_wp ! - -69 70 68 !!---------------------------------------------------------------------- 71 69 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 92 90 ! 93 91 INTEGER :: ji, jj, jk, jl ! dummy loop indices 94 REAL(wp) :: zinda, zindb95 92 !!------------------------------------------------------------------ 96 93 … … 111 108 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 112 109 ! 113 zinda= MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )114 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda! ice thickness110 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 111 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch ! ice thickness 115 112 END DO 116 113 END DO … … 132 129 DO jj = 1, jpj 133 130 DO ji = 1, jpi 134 zinda = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )135 zindb = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )136 131 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 137 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * zinda ! ice salinity 138 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) * zindb ! ice age 132 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 133 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch ! ice salinity 134 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 135 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) * rswitch ! ice age 139 136 END DO 140 137 END DO … … 161 158 INTEGER :: ji, jj, jk, jl ! dummy loop indices 162 159 REAL(wp) :: zq_i, zaaa, zbbb, zccc, zdiscrim ! local scalars 163 REAL(wp) :: ztmelts, z indb, zq_s, zfac1, zfac2 ! - -160 REAL(wp) :: ztmelts, zq_s, zfac1, zfac2 ! - - 164 161 !!------------------------------------------------------------------ 165 162 … … 170 167 DO jj = 1, jpj 171 168 DO ji = 1, jpi 172 zindb= 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes173 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb174 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb175 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb169 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 170 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 171 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 172 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 176 173 END DO 177 174 END DO … … 182 179 DO jj = 1, jpj 183 180 DO ji = 1, jpi 184 zindb= 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes185 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * zindb181 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 182 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * rswitch 186 183 END DO 187 184 END DO … … 203 200 DO ji = 1, jpi 204 201 ! ! Energy of melting q(S,T) [J.m-3] 205 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) ! zindb= 0 if no ice and 1 if yes206 zq_i = zindb* e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp)202 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) ! rswitch = 0 if no ice and 1 if yes 203 zq_i = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 207 204 zq_i = zq_i * unit_fac !convert units 208 205 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt ! Ice layer melt temperature … … 212 209 zccc = lfus * (ztmelts-rtt) 213 210 zdiscrim = SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 214 t_i(ji,jj,jk,jl) = rtt + zindb*( - zbbb - zdiscrim ) / ( 2.0 *zaaa )211 t_i(ji,jj,jk,jl) = rtt + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) 215 212 t_i(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) ) ! 100-rtt < t_i < rtt 216 213 END DO … … 229 226 DO ji = 1, jpi 230 227 !Energy of melting q(S,T) [J.m-3] 231 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) ) ! zindb= 0 if no ice and 1 if yes232 zq_s = zindb* e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp)228 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) ) ! rswitch = 0 if no ice and 1 if yes 229 zq_s = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp) 233 230 zq_s = zq_s * unit_fac ! convert units 234 231 ! 235 t_s(ji,jj,jk,jl) = rtt + zindb* ( - zfac1 * zq_s + zfac2 )232 t_s(ji,jj,jk,jl) = rtt + rswitch * ( - zfac1 * zq_s + zfac2 ) 236 233 t_s(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15, t_s(ji,jj,jk,jl) ) ) ! 100-rtt < t_i < rtt 237 234 END DO … … 248 245 DO jj = 1, jpj 249 246 DO ji = 1, jpi 250 zindb= ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) )251 tm_i(ji,jj) = tm_i(ji,jj) + zindb* t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) &247 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 248 tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 252 249 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 253 250 END DO … … 295 292 INTEGER :: ji, jj, jk, jl ! dummy loop index 296 293 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac, zsal ! local scalar 297 REAL(wp) :: z ind0, zind01, zindbal, zargtemp , zs_zero ! - -294 REAL(wp) :: zswi0, zswi01, zswibal, zargtemp , zs_zero ! - - 298 295 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_slope_s, zalpha ! 3D pointer 299 296 !!------------------------------------------------------------------ … … 330 327 DO jj = 1, jpj 331 328 DO ji = 1, jpi 332 ! z ind0 = 1 if sm_i le s_i_0 and 0 otherwise333 z ind0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i(ji,jj,jl) ) )334 ! z ind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws335 z ind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i(ji,jj,jl) ) )336 ! If 2.sm_i GE sss_m then z indbal = 1329 ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise 330 zswi0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i(ji,jj,jl) ) ) 331 ! zswi01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws 332 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i(ji,jj,jl) ) ) 333 ! If 2.sm_i GE sss_m then zswibal = 1 337 334 ! this is to force a constant salinity profile in the Baltic Sea 338 z indbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )339 zalpha(ji,jj,jl) = z ind0 + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )340 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - z indbal )335 zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 336 zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 ) 337 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zswibal ) 341 338 END DO 342 339 END DO … … 390 387 !!------------------------------------------------------------------ 391 388 INTEGER :: ji, jj, jk, jl ! dummy loop indices 392 REAL(wp) :: zindb ! - -393 389 !!------------------------------------------------------------------ 394 390 … … 399 395 DO jj = 1, jpj 400 396 DO ji = 1, jpi 401 zindb= ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) )402 tm_i(ji,jj) = tm_i(ji,jj) + zindb* t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) &397 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 398 tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 403 399 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 404 400 END DO … … 421 417 !!------------------------------------------------------------------ 422 418 INTEGER :: ji, jj, jk, jl ! dummy loop indices 423 REAL(wp) :: zbvi , zinda, zindb! local scalars419 REAL(wp) :: zbvi ! local scalars 424 420 !!------------------------------------------------------------------ 425 421 ! … … 429 425 DO jj = 1, jpj 430 426 DO ji = 1, jpi 431 zinda = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) ) ) 432 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 433 zbvi = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 ) & 427 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) ) ) 428 zbvi = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 ) & 434 429 & * v_i(ji,jj,jl) / REAL(nlay_i,wp) 435 bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi / MAX( vt_i(ji,jj) , epsi10 ) 430 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 431 bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi / MAX( vt_i(ji,jj) , epsi10 ) 436 432 END DO 437 433 END DO … … 454 450 INTEGER :: ii, ij ! local integers 455 451 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal ! local scalars 456 REAL(wp) :: zalpha, z ind0, zind01, zindbal, zs_zero ! - -452 REAL(wp) :: zalpha, zswi0, zswi01, zswibal, zs_zero ! - - 457 453 ! 458 454 REAL(wp), POINTER, DIMENSION(:) :: z_slope_s … … 488 484 ii = MOD( npb(ji) - 1 , jpi ) + 1 489 485 ij = ( npb(ji) - 1 ) / jpi + 1 490 ! z ind0 = 1 if sm_i le s_i_0 and 0 otherwise491 z ind0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i_1d(ji) ) )492 ! z ind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws493 z ind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_1d(ji) ) )494 ! if 2.sm_i GE sss_m then z indbal = 1486 ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise 487 zswi0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i_1d(ji) ) ) 488 ! zswi01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws 489 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_1d(ji) ) ) 490 ! if 2.sm_i GE sss_m then zswibal = 1 495 491 ! this is to force a constant salinity profile in the Baltic Sea 496 z indbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )492 zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 497 493 ! 498 zalpha = ( z ind0 + zind01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zindbal )494 zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zswibal ) 499 495 ! 500 496 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2 -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r4967 r4972 35 35 PUBLIC lim_wri_state ! called by dia_wri_state 36 36 37 REAL(wp) :: epsi06 = 1.e-6_wp38 37 !!---------------------------------------------------------------------- 39 38 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 59 58 INTEGER, INTENT(in) :: kindic ! if kindic < 0 there has been an error somewhere 60 59 ! 61 INTEGER :: ji, jj, jk, jl ! dummy loop indices62 REAL(wp) :: z inda, zindb, z1_36560 INTEGER :: ji, jj, jk, jl ! dummy loop indices 61 REAL(wp) :: z1_365 63 62 REAL(wp) :: ztmp 64 REAL(wp), POINTER, DIMENSION(:,:,:) :: 65 REAL(wp), POINTER, DIMENSION(:,:) :: z2d, z2da, z2db, zind! 2D workspace63 REAL(wp), POINTER, DIMENSION(:,:,:) :: zoi, zei 64 REAL(wp), POINTER, DIMENSION(:,:) :: z2d, z2da, z2db, zswi ! 2D workspace 66 65 !!------------------------------------------------------------------- 67 66 … … 69 68 70 69 CALL wrk_alloc( jpi, jpj, jpl, zoi, zei ) 71 CALL wrk_alloc( jpi, jpj , z2d, z2da, z2db, z ind)70 CALL wrk_alloc( jpi, jpj , z2d, z2da, z2db, zswi ) 72 71 73 72 !----------------------------- … … 81 80 DO jj = 1, jpj ! presence indicator of ice 82 81 DO ji = 1, jpi 83 z ind(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) )82 zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) 84 83 END DO 85 84 END DO … … 90 89 DO jj = 1, jpj 91 90 DO ji = 1, jpi 92 z2d(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * z ind(ji,jj)91 z2d(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) 93 92 END DO 94 93 END DO … … 99 98 DO jj = 1, jpj 100 99 DO ji = 1, jpi 101 z2d(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * z ind(ji,jj)100 z2d(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) 102 101 END DO 103 102 END DO … … 129 128 DO jj = 1, jpj 130 129 DO ji = 1, jpi 131 z2d(ji,jj) = z2d(ji,jj) + z ind(ji,jj) * oa_i(ji,jj,jl)130 z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * oa_i(ji,jj,jl) 132 131 END DO 133 132 END DO … … 140 139 DO jj = 1, jpj 141 140 DO ji = 1, jpi 142 z2d(ji,jj) = ( tm_i(ji,jj) - rtt ) * z ind(ji,jj)141 z2d(ji,jj) = ( tm_i(ji,jj) - rtt ) * zswi(ji,jj) 143 142 END DO 144 143 END DO … … 151 150 DO jj = 1, jpj 152 151 DO ji = 1, jpi 153 z2d(ji,jj) = z2d(ji,jj) + z ind(ji,jj) * ( t_su(ji,jj,jl) - rtt ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 )152 z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * ( t_su(ji,jj,jl) - rtt ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 ) 154 153 END DO 155 154 END DO … … 161 160 DO jj = 1, jpj 162 161 DO ji = 1, jpi 163 zindb= MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) ) )164 z2d(ji,jj) = hicol(ji,jj) * zindb162 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) ) ) 163 z2d(ji,jj) = hicol(ji,jj) * rswitch 165 164 END DO 166 165 END DO … … 245 244 DO jj = 1, jpj 246 245 DO ji = 1, jpi 247 zinda= MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) )248 zoi(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi06 ) * zinda246 rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 247 zoi(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi06 ) * rswitch 249 248 END DO 250 249 END DO … … 260 259 DO jj = 1, jpj 261 260 DO ji = 1, jpi 262 zinda= MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) )261 rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 263 262 zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0* & 264 263 ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rtt ), - epsi06 ) ) * & 265 zinda/ nlay_i264 rswitch / nlay_i 266 265 END DO 267 266 END DO … … 276 275 277 276 CALL wrk_dealloc( jpi, jpj, jpl, zoi, zei ) 278 CALL wrk_dealloc( jpi, jpj , z2d, z ind, z2da, z2db )277 CALL wrk_dealloc( jpi, jpj , z2d, zswi, z2da, z2db ) 279 278 280 279 IF( nn_timing == 1 ) CALL timing_stop('limwri') -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90
r4924 r4972 45 45 PUBLIC bdy_ice_lim_dyn ! routine called in limrhg 46 46 47 REAL(wp) :: epsi20 = 1.e-20_wp ! module constants48 REAL(wp) :: epsi10 = 1.e-10_wp ! min area allowed by ice model49 47 !!---------------------------------------------------------------------- 50 48 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 98 96 INTEGER :: ji, jj, ii, ij ! local scalar 99 97 REAL(wp) :: zwgt, zwgt1 ! local scalar 100 REAL(wp) :: z inda, ztmelts, zdh98 REAL(wp) :: ztmelts, zdh 101 99 102 100 !!------------------------------------------------------------------------------ … … 115 113 hicif(ji,jj) = ( hicif(ji,jj) * zwgt1 + dta%hicif(jb) * zwgt ) * tmask(ji,jj,1) ! Ice depth 116 114 hsnif(ji,jj) = ( hsnif(ji,jj) * zwgt1 + dta%hsnif(jb) * zwgt ) * tmask(ji,jj,1) ! Snow depth 117 118 ! zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - frld(ji,jj) ) ) ! 0 if no ice119 ! !------------------------------120 ! ! Sea ice surface temperature121 ! !------------------------------122 ! sist(ji,jj) = zinda * 270.0 + ( 1.0 - zinda ) * tfu(ji,jj)123 ! !-----------------------------------------------124 ! ! Ice/snow temperatures and energy stored in brines125 ! !-----------------------------------------------126 ! !!! TO BE CONTIUNED (as LIM3 below) !!!127 ! zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) ) ! = 0 for SH, =1 for NH128 !129 ! ! Recover in situ values.130 ! zindb = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) )131 ! zacrith = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) )132 ! zs0a (ji,jj) = zindb * MIN( zs0a(ji,jj), zacrith )133 ! hsnif(ji,jj) = zindb * ( zs0sn(ji,jj) /MAX( zs0a(ji,jj), epsi16 ) )134 ! hicif(ji,jj) = zindb * ( zs0ice(ji,jj)/MAX( zs0a(ji,jj), epsi16 ) )135 ! zindsn = MAX( rzero, SIGN( rone, hsnif(ji,jj) - epsi06 ) )136 ! zindic = MAX( rzero, SIGN( rone, hicif(ji,jj) - epsi03 ) )137 ! zindb = MAX( zindsn, zindic )138 ! zs0a (ji,jj) = zindb * zs0a(ji,jj)139 ! frld (ji,jj) = 1.0 - zs0a(ji,jj)140 ! hsnif(ji,jj) = zindsn * hsnif(ji,jj)141 ! hicif(ji,jj) = zindic * hicif(ji,jj)142 ! zusvosn = 1.0/MAX( hsnif(ji,jj) * zs0a(ji,jj), epsi16 )143 ! zusvoic = 1.0/MAX( hicif(ji,jj) * zs0a(ji,jj), epsi16 )144 ! zignm = MAX( rzero, SIGN( rone, hsndif - hsnif(ji,jj) ) )145 ! zrtt = 173.15 * rone146 ! ztsn = zignm * tbif(ji,jj,1) &147 ! + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) )148 ! ztic1 = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) )149 ! ztic2 = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) )150 !151 ! tbif(ji,jj,1) = zindsn * ztsn + ( 1.0 - zindsn ) * tfu(ji,jj)152 ! tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj)153 ! tbif(ji,jj,3) = zindic * ztic2 + ( 1.0 - zindic ) * tfu(ji,jj)154 ! qstoif(ji,jj) = zindb * xlic * zs0st(ji,jj) / MAX( zs0a(ji,jj), epsi16 )155 156 115 END DO 157 116 … … 215 174 IF ( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj+1,1) == 0. ) jpbound = 1; ii = ji ; ij = jj-1 216 175 217 zinda= 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ii,ij) + 0.01 ) ) ! 0 if no ice176 rswitch = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ii,ij) + 0.01 ) ) ! 0 if no ice 218 177 219 178 ! concentration and thickness 220 a_i (ji,jj,jl) = a_i (ii,ij,jl) * zinda221 ht_i(ji,jj,jl) = ht_i(ii,ij,jl) * zinda222 ht_s(ji,jj,jl) = ht_s(ii,ij,jl) * zinda179 a_i (ji,jj,jl) = a_i (ii,ij,jl) * rswitch 180 ht_i(ji,jj,jl) = ht_i(ii,ij,jl) * rswitch 181 ht_s(ji,jj,jl) = ht_s(ii,ij,jl) * rswitch 223 182 224 183 ! Ice and snow volumes … … 231 190 232 191 ! Ice salinity, age, temperature 233 sm_i(ji,jj,jl) = zinda * rn_ice_sal(ib_bdy) + ( 1.0 - zinda) * s_i_min234 o_i(ji,jj,jl) = zinda * rn_ice_age(ib_bdy) + ( 1.0 - zinda)235 t_su(ji,jj,jl) = zinda * rn_ice_tem(ib_bdy) + ( 1.0 - zinda) * rn_ice_tem(ib_bdy)192 sm_i(ji,jj,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * s_i_min 193 o_i(ji,jj,jl) = rswitch * rn_ice_age(ib_bdy) + ( 1.0 - rswitch ) 194 t_su(ji,jj,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rn_ice_tem(ib_bdy) 236 195 DO jk = 1, nlay_s 237 t_s(ji,jj,jk,jl) = zinda * rn_ice_tem(ib_bdy) + ( 1.0 - zinda) * rtt196 t_s(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rtt 238 197 END DO 239 198 DO jk = 1, nlay_i 240 t_i(ji,jj,jk,jl) = zinda * rn_ice_tem(ib_bdy) + ( 1.0 - zinda) * rtt241 s_i(ji,jj,jk,jl) = zinda * rn_ice_sal(ib_bdy) + ( 1.0 - zinda) * s_i_min199 t_i(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rtt 200 s_i(ji,jj,jk,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * s_i_min 242 201 END DO 243 202 … … 245 204 246 205 ! Ice salinity, age, temperature 247 sm_i(ji,jj,jl) = zinda * sm_i(ii,ij,jl) + ( 1.0 - zinda) * s_i_min248 o_i(ji,jj,jl) = zinda * o_i(ii,ij,jl) + ( 1.0 - zinda)249 t_su(ji,jj,jl) = zinda * t_su(ii,ij,jl) + ( 1.0 - zinda) * rtt206 sm_i(ji,jj,jl) = rswitch * sm_i(ii,ij,jl) + ( 1.0 - rswitch ) * s_i_min 207 o_i(ji,jj,jl) = rswitch * o_i(ii,ij,jl) + ( 1.0 - rswitch ) 208 t_su(ji,jj,jl) = rswitch * t_su(ii,ij,jl) + ( 1.0 - rswitch ) * rtt 250 209 DO jk = 1, nlay_s 251 t_s(ji,jj,jk,jl) = zinda * t_s(ii,ij,jk,jl) + ( 1.0 - zinda) * rtt210 t_s(ji,jj,jk,jl) = rswitch * t_s(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rtt 252 211 END DO 253 212 DO jk = 1, nlay_i 254 t_i(ji,jj,jk,jl) = zinda * t_i(ii,ij,jk,jl) + ( 1.0 - zinda) * rtt255 s_i(ji,jj,jk,jl) = zinda * s_i(ii,ij,jk,jl) + ( 1.0 - zinda) * s_i_min213 t_i(ji,jj,jk,jl) = rswitch * t_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rtt 214 s_i(ji,jj,jk,jl) = rswitch * s_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * s_i_min 256 215 END DO 257 216 … … 269 228 DO jk = 1, nlay_s 270 229 ! Snow energy of melting 271 e_s(ji,jj,jk,jl) = zinda* rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )230 e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 272 231 ! Change dimensions 273 232 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac … … 278 237 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 279 238 ! heat content per unit volume 280 e_i(ji,jj,jk,jl) = zinda* rhoic * &239 e_i(ji,jj,jk,jl) = rswitch * rhoic * & 281 240 ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 282 241 + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & … … 335 294 INTEGER :: ji, jj ! local scalar 336 295 INTEGER :: ib_bdy ! Loop index 337 REAL(wp) :: zmsk1, zmsk2, zflag , zinda296 REAL(wp) :: zmsk1, zmsk2, zflag 338 297 !!------------------------------------------------------------------------------ 339 298 ! … … 375 334 ENDIF 376 335 ! mask ice velocities 377 zinda= 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice378 u_ice(ji,jj) = zinda* u_ice(ji,jj)336 rswitch = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice 337 u_ice(ji,jj) = rswitch * u_ice(ji,jj) 379 338 380 339 ENDDO … … 404 363 ENDIF 405 364 ! mask ice velocities 406 zinda= 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice407 v_ice(ji,jj) = zinda* v_ice(ji,jj)365 rswitch = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice 366 v_ice(ji,jj) = rswitch * v_ice(ji,jj) 408 367 409 368 ENDDO -
branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r4946 r4972 221 221 v_ice_b(:,:) = v_ice(:,:) 222 222 223 ! trends !!gm is it truly necessary ???224 d_a_i_thd (:,:,:) = 0._wp ; d_a_i_trp (:,:,:) = 0._wp225 d_v_i_thd (:,:,:) = 0._wp ; d_v_i_trp (:,:,:) = 0._wp226 d_e_i_thd (:,:,:,:) = 0._wp ; d_e_i_trp (:,:,:,:) = 0._wp227 d_v_s_thd (:,:,:) = 0._wp ; d_v_s_trp (:,:,:) = 0._wp228 d_e_s_thd (:,:,:,:) = 0._wp ; d_e_s_trp (:,:,:,:) = 0._wp229 d_smv_i_thd(:,:,:) = 0._wp ; d_smv_i_trp(:,:,:) = 0._wp230 d_oa_i_thd (:,:,:) = 0._wp ; d_oa_i_trp (:,:,:) = 0._wp231 d_u_ice_dyn(:,:) = 0._wp ; d_v_ice_dyn(:,:) = 0._wp232 233 223 ! salt, heat and mass fluxes 234 224 sfx (:,:) = 0._wp ; … … 254 244 hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp 255 245 hfx_err(:,:) = 0._wp ; hfx_err_rem(:,:) = 0._wp 256 257 !258 fhld (:,:) = 0._wp259 fmmflx(:,:) = 0._wp260 ! part of solar radiation transmitted through the ice261 ftr_ice(:,:,:) = 0._wp262 263 ! diags264 diag_trp_vi (:,:) = 0._wp ; diag_trp_vs(:,:) = 0._wp ; diag_trp_ei(:,:) = 0._wp ; diag_trp_es(:,:) = 0._wp265 diag_heat_dhc(:,:) = 0._wp266 267 ! dynamical invariants268 delta_i(:,:) = 0._wp ; divu_i(:,:) = 0._wp ; shear_i(:,:) = 0._wp269 246 270 247 CALL lim_rst_opn( kt ) ! Open Ice restart file
Note: See TracChangeset
for help on using the changeset viewer.
