Changeset 5064 for branches/2015/dev_r5044_CNRS_LIM3CLEAN
- Timestamp:
- 2015-02-05T18:54:24+01:00 (9 years ago)
- Location:
- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO
- Files:
-
- 22 edited
-
LIM_SRC_3/dom_ice.F90 (modified) (2 diffs)
-
LIM_SRC_3/ice.F90 (modified) (3 diffs)
-
LIM_SRC_3/limadv.F90 (modified) (2 diffs)
-
LIM_SRC_3/limcons.F90 (modified) (3 diffs)
-
LIM_SRC_3/limctl.F90 (modified) (4 diffs)
-
LIM_SRC_3/limdiahsb.F90 (modified) (1 diff)
-
LIM_SRC_3/limdyn.F90 (modified) (2 diffs)
-
LIM_SRC_3/limistate.F90 (modified) (2 diffs)
-
LIM_SRC_3/limitd_me.F90 (modified) (4 diffs)
-
LIM_SRC_3/limitd_th.F90 (modified) (2 diffs)
-
LIM_SRC_3/limmsh.F90 (modified) (1 diff)
-
LIM_SRC_3/limrhg.F90 (modified) (18 diffs)
-
LIM_SRC_3/limsbc.F90 (modified) (3 diffs)
-
LIM_SRC_3/limthd.F90 (modified) (12 diffs)
-
LIM_SRC_3/limthd_dif.F90 (modified) (3 diffs)
-
LIM_SRC_3/limthd_ent.F90 (modified) (1 diff)
-
LIM_SRC_3/limthd_lac.F90 (modified) (2 diffs)
-
LIM_SRC_3/limtrp.F90 (modified) (5 diffs)
-
LIM_SRC_3/limupdate1.F90 (modified) (4 diffs)
-
LIM_SRC_3/limupdate2.F90 (modified) (4 diffs)
-
LIM_SRC_3/limvar.F90 (modified) (4 diffs)
-
OPA_SRC/BDY/bdyice_lim.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90
r5059 r5064 21 21 22 22 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fcor !: coriolis coefficient 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: area !: surface of grid cell24 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tms, tmi !: temperature mask, mask for stress 25 24 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmu, tmv !: mask at u and v velocity points … … 42 41 !!------------------------------------------------------------------- 43 42 ! 44 ALLOCATE( fcor(jpi,jpj) , area(jpi,jpj) ,&43 ALLOCATE( fcor(jpi,jpj) , & 45 44 & tms (jpi,jpj) , tmi (jpi,jpj) , & 46 45 & tmu (jpi,jpj) , tmv (jpi,jpj) , & -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r5059 r5064 109 109 !! smv_i | - | Sea ice salt content | ppt.m | 110 110 !! oa_i ! - ! Sea ice areal age content | day | 111 !! e_i ! - ! Ice enthalpy | 10^9 J|111 !! e_i ! - ! Ice enthalpy | J/m2 | 112 112 !! - ! q_i_1d ! Ice enthalpy per unit vol. | J/m3 | 113 !! e_s ! - ! Snow enthalpy | 10^9 J|113 !! e_s ! - ! Snow enthalpy | J/m2 | 114 114 !! - ! q_s_1d ! Snow enthalpy per unit vol. | J/m3 | 115 115 !! | … … 224 224 225 225 ! !!** define some parameters 226 REAL(wp), PUBLIC, PARAMETER :: unit_fac = 1.e+09_wp !: conversion factor for ice / snow enthalpy227 226 REAL(wp), PUBLIC, PARAMETER :: epsi06 = 1.e-06_wp !: small number 228 227 REAL(wp), PUBLIC, PARAMETER :: epsi10 = 1.e-10_wp !: small number … … 331 330 332 331 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_i !: ice temperatures [K] 333 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i !: ice thermal contents [Giga J]332 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i !: ice thermal contents [J/m2] 334 333 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: s_i !: ice salinities [PSU] 335 334 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
r4990 r5064 97 97 98 98 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 99 psm (:,:) = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )99 psm (:,:) = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 100 100 101 101 ! Calculate fluxes and moments between boxes i<-->i+1 … … 282 282 283 283 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 284 psm(:,:) = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )284 psm(:,:) = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 285 285 286 286 ! Calculate fluxes and moments between boxes j<-->j+1 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r5059 r5064 167 167 REAL(wp) :: zvi, zsmv, zei, zfs, zfw, zft 168 168 REAL(wp) :: zvmin, zamin, zamax 169 REAL(wp) :: zconv 170 171 zconv = 1.e-9 169 172 170 173 IF( icount == 0 ) THEN 171 174 172 zvi_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 173 zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 174 zei_b = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 175 zvi_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tms(:,:) ) 176 zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tms(:,:) ) 177 zei_b = glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv + & 178 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv ) 175 179 zfw_b = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 176 180 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 177 & ) * area(:,:) * tms(:,:) )181 & ) * e12t(:,:) * tms(:,:) ) 178 182 zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 179 183 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 180 & ) * area(:,:) * tms(:,:) )184 & ) * e12t(:,:) * tms(:,:) ) 181 185 zft_b = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 182 186 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 183 & ) * area(:,:) / unit_fac * tms(:,:))187 & ) * e12t(:,:) * tms(:,:) * zconv ) 184 188 185 189 ELSEIF( icount == 1 ) THEN … … 187 191 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 188 192 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 189 & ) * area(:,:) * tms(:,:) ) - zfs_b193 & ) * e12t(:,:) * tms(:,:) ) - zfs_b 190 194 zfw = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 191 195 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 192 & ) * area(:,:) * tms(:,:) ) - zfw_b196 & ) * e12t(:,:) * tms(:,:) ) - zfw_b 193 197 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 194 198 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 195 & ) * area(:,:) / unit_fac * tms(:,:)) - zft_b199 & ) * e12t(:,:) * tms(:,:) * zconv ) - zft_b 196 200 197 zvi = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw 198 zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 199 zei = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zei_b * r1_rdtice + zft 201 zvi = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw 202 zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 203 zei = glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv + & 204 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft 200 205 201 206 zvmin = glob_min(v_i) … … 206 211 IF ( ABS( zvi ) > 1.e-4 ) WRITE(numout,*) 'violation volume [kg/day] (',cd_routine,') = ',(zvi * rday) 207 212 IF ( ABS( zsmv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday) 208 IF ( ABS( zei ) > 1. ) WRITE(numout,*) 'violation enthalpy [1e9 J](',cd_routine,') = ',(zei)209 IF ( zvmin < 1.e-10) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',(zvmin)210 IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+ 1.e-10 ) THEN213 IF ( ABS( zei ) > 1.e-2 ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',(zei) 214 IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',(zvmin) 215 IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+epsi10 ) THEN 211 216 WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax 212 217 ENDIF 213 IF ( zamin < 1.e-10) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin218 IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin 214 219 ENDIF 215 220 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
r5059 r5064 306 306 WRITE(numout,*) ' - Cell values ' 307 307 WRITE(numout,*) ' ~~~~~~~~~~~ ' 308 WRITE(numout,*) ' cell area : ', area(ji,jj)308 WRITE(numout,*) ' cell area : ', e12t(ji,jj) 309 309 WRITE(numout,*) ' at_i : ', at_i(ji,jj) 310 310 WRITE(numout,*) ' vt_i : ', vt_i(ji,jj) … … 317 317 WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl) 318 318 WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl) 319 WRITE(numout,*) ' e_s : ', e_s(ji,jj,1,jl) /1.0e9320 WRITE(numout,*) ' e_i : ', e_i(ji,jj,1:nlay_i,jl) /1.0e9319 WRITE(numout,*) ' e_s : ', e_s(ji,jj,1,jl) 320 WRITE(numout,*) ' e_i : ', e_i(ji,jj,1:nlay_i,jl) 321 321 WRITE(numout,*) ' t_su : ', t_su(ji,jj,jl) 322 322 WRITE(numout,*) ' t_snow : ', t_s(ji,jj,1,jl) … … 350 350 WRITE(numout,*) ' - Cell values ' 351 351 WRITE(numout,*) ' ~~~~~~~~~~~ ' 352 WRITE(numout,*) ' cell area : ', area(ji,jj)352 WRITE(numout,*) ' cell area : ', e12t(ji,jj) 353 353 WRITE(numout,*) ' at_i : ', at_i(ji,jj) 354 354 WRITE(numout,*) ' vt_i : ', vt_i(ji,jj) … … 372 372 WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl) , ' v_i_b : ', v_i_b(ji,jj,jl) 373 373 WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl) , ' v_s_b : ', v_s_b(ji,jj,jl) 374 WRITE(numout,*) ' e_i1 : ', e_i(ji,jj,1,jl) /1.0e9 , ' ei1 : ', e_i_b(ji,jj,1,jl)/1.0e9375 WRITE(numout,*) ' e_i2 : ', e_i(ji,jj,2,jl) /1.0e9 , ' ei2_b : ', e_i_b(ji,jj,2,jl)/1.0e9374 WRITE(numout,*) ' e_i1 : ', e_i(ji,jj,1,jl) , ' ei1 : ', e_i_b(ji,jj,1,jl) 375 WRITE(numout,*) ' e_i2 : ', e_i(ji,jj,2,jl) , ' ei2_b : ', e_i_b(ji,jj,2,jl) 376 376 WRITE(numout,*) ' e_snow : ', e_s(ji,jj,1,jl) , ' e_snow_b : ', e_s_b(ji,jj,1,jl) 377 377 WRITE(numout,*) ' smv_i : ', smv_i(ji,jj,jl) , ' smv_i_b : ', smv_i_b(ji,jj,jl) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
r5055 r5064 71 71 72 72 ! 1/area 73 z1_area = 1._wp / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 )74 75 rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) )73 z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tms(:,:) ), epsi06 ) 74 75 rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tms(:,:) ) - epsi06 ) ) 76 76 ! ----------------------- ! 77 77 ! 1 - Content variations ! 78 78 ! ----------------------- ! 79 zbg_ivo = glob_sum( vt_i(:,:) * area(:,:) * tms(:,:) ) ! volume ice80 zbg_svo = glob_sum( vt_s(:,:) * area(:,:) * tms(:,:) ) ! volume snow81 zbg_are = glob_sum( at_i(:,:) * area(:,:) * tms(:,:) ) ! area82 zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) ! mean salt content83 zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * area(:,:) * tms(:,:) ) ! mean temp content84 85 !zbg_ihc = glob_sum( et_i(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content86 !zbg_shc = glob_sum( et_s(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content79 zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tms(:,:) ) ! volume ice 80 zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tms(:,:) ) ! volume snow 81 zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tms(:,:) ) ! area 82 zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tms(:,:) ) ! mean salt content 83 zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * e12t(:,:) * tms(:,:) ) ! mean temp content 84 85 !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tms(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content 86 !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tms(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content 87 87 88 88 ! Volume 89 89 ztmp = rswitch * z1_area * r1_rau0 * rday 90 zbg_vfx = ztmp * glob_sum( emp(:,:) * area(:,:) * tms(:,:) )91 zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * area(:,:) * tms(:,:) )92 zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * area(:,:) * tms(:,:) )93 zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * area(:,:) * tms(:,:) )94 zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * area(:,:) * tms(:,:) )95 zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * area(:,:) * tms(:,:) )96 zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * area(:,:) * tms(:,:) )97 zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * area(:,:) * tms(:,:) )98 zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * area(:,:) * tms(:,:) )99 zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * area(:,:) * tms(:,:) )100 zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * area(:,:) * tms(:,:) )90 zbg_vfx = ztmp * glob_sum( emp(:,:) * e12t(:,:) * tms(:,:) ) 91 zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tms(:,:) ) 92 zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tms(:,:) ) 93 zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tms(:,:) ) 94 zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tms(:,:) ) 95 zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tms(:,:) ) 96 zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tms(:,:) ) 97 zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tms(:,:) ) 98 zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tms(:,:) ) 99 zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tms(:,:) ) 100 zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tms(:,:) ) 101 101 102 102 ! Salt 103 zbg_sfx = ztmp * glob_sum( sfx(:,:) * area(:,:) * tms(:,:) )104 zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) )105 zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) )106 zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * area(:,:) * tms(:,:) )107 108 zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * area(:,:) * tms(:,:) )109 zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * area(:,:) * tms(:,:) )110 zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * area(:,:) * tms(:,:) )111 zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * area(:,:) * tms(:,:) )112 zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * area(:,:) * tms(:,:) )103 zbg_sfx = ztmp * glob_sum( sfx(:,:) * e12t(:,:) * tms(:,:) ) 104 zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tms(:,:) ) 105 zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tms(:,:) ) 106 zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tms(:,:) ) 107 108 zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tms(:,:) ) 109 zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tms(:,:) ) 110 zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tms(:,:) ) 111 zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tms(:,:) ) 112 zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tms(:,:) ) 113 113 114 114 ! Heat budget 115 115 zbg_ihc = glob_sum( et_i(:,:) * 1.e-20 ) ! ice heat content [1.e-20 J] 116 116 zbg_shc = glob_sum( et_s(:,:) * 1.e-20 ) ! snow heat content [1.e-20 J] 117 zbg_hfx_dhc = glob_sum( diag_heat_dhc(:,:) * area(:,:) * tms(:,:) ) ! [in W]118 zbg_hfx_spr = glob_sum( hfx_spr(:,:) * area(:,:) * tms(:,:) ) ! [in W]119 120 zbg_hfx_thd = glob_sum( hfx_thd(:,:) * area(:,:) * tms(:,:) ) ! [in W]121 zbg_hfx_dyn = glob_sum( hfx_dyn(:,:) * area(:,:) * tms(:,:) ) ! [in W]122 zbg_hfx_res = glob_sum( hfx_res(:,:) * area(:,:) * tms(:,:) ) ! [in W]123 zbg_hfx_sub = glob_sum( hfx_sub(:,:) * area(:,:) * tms(:,:) ) ! [in W]124 zbg_hfx_snw = glob_sum( hfx_snw(:,:) * area(:,:) * tms(:,:) ) ! [in W]125 zbg_hfx_sum = glob_sum( hfx_sum(:,:) * area(:,:) * tms(:,:) ) ! [in W]126 zbg_hfx_bom = glob_sum( hfx_bom(:,:) * area(:,:) * tms(:,:) ) ! [in W]127 zbg_hfx_bog = glob_sum( hfx_bog(:,:) * area(:,:) * tms(:,:) ) ! [in W]128 zbg_hfx_dif = glob_sum( hfx_dif(:,:) * area(:,:) * tms(:,:) ) ! [in W]129 zbg_hfx_opw = glob_sum( hfx_opw(:,:) * area(:,:) * tms(:,:) ) ! [in W]130 zbg_hfx_out = glob_sum( hfx_out(:,:) * area(:,:) * tms(:,:) ) ! [in W]131 zbg_hfx_in = glob_sum( hfx_in(:,:) * area(:,:) * tms(:,:) ) ! [in W]117 zbg_hfx_dhc = glob_sum( diag_heat_dhc(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 118 zbg_hfx_spr = glob_sum( hfx_spr(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 119 120 zbg_hfx_thd = glob_sum( hfx_thd(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 121 zbg_hfx_dyn = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 122 zbg_hfx_res = glob_sum( hfx_res(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 123 zbg_hfx_sub = glob_sum( hfx_sub(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 124 zbg_hfx_snw = glob_sum( hfx_snw(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 125 zbg_hfx_sum = glob_sum( hfx_sum(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 126 zbg_hfx_bom = glob_sum( hfx_bom(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 127 zbg_hfx_bog = glob_sum( hfx_bog(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 128 zbg_hfx_dif = glob_sum( hfx_dif(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 129 zbg_hfx_opw = glob_sum( hfx_opw(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 130 zbg_hfx_out = glob_sum( hfx_out(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 131 zbg_hfx_in = glob_sum( hfx_in(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 132 132 133 133 ! --------------------------------------------- ! 134 134 ! 2 - Trends due to forcing and ice growth/melt ! 135 135 ! --------------------------------------------- ! 136 z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * area(:,:) * tms(:,:) ) ! volume fluxes137 z_frc_sal = r1_rau0 * glob_sum( sfx(:,:) * area(:,:) * tms(:,:) ) ! salt fluxes136 z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tms(:,:) ) ! volume fluxes 137 z_frc_sal = r1_rau0 * glob_sum( sfx(:,:) * e12t(:,:) * tms(:,:) ) ! salt fluxes 138 138 z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 139 & wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * area(:,:) * tms(:,:) ) ! volume fluxes139 & wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * e12t(:,:) * tms(:,:) ) ! volume fluxes 140 140 ! 141 141 frc_vol = frc_vol + z_frc_vol * rdt_ice -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r5059 r5064 191 191 CALL prt_ctl(tab2d_1=delta_i , clinfo1=' lim_dyn : delta_i :') 192 192 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_dyn : strength :') 193 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_dyn : cell area :')193 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_dyn : cell area :') 194 194 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_dyn : at_i :') 195 195 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_dyn : vt_i :') … … 274 274 rhoco = rau0 * cw 275 275 ! 276 ! Diffusion coefficients .276 ! Diffusion coefficients 277 277 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 278 278 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 279 279 280 280 za00 = ahi0 / ( 1.e05_wp ) ! 1.e05 = 100km = max grid space at 60° latitude in orca2 281 ! (i.e.60° = min latitude for ice cover)281 ! (60° = min latitude for ice cover) 282 282 DO jj = 1, jpj 283 283 DO ji = 1, jpi -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r5055 r5064 345 345 ! Snow energy of melting 346 346 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 347 ! Change dimensions 348 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 349 ! Multiply by volume, so that heat content in Joules 350 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 347 348 ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 349 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) / nlay_s 351 350 END DO ! ji 352 351 END DO ! jj … … 368 367 - rcp * ( ztmelts - rtt ) ) 369 368 370 ! Correct dimensions to avoid big values 371 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 372 373 ! Mutliply by ice volume, and divide by number of layers to get heat content in J 374 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i 369 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 370 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / nlay_i 375 371 END DO ! ji 376 372 END DO ! jj -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r5055 r5064 339 339 !-----------------------------------------------------------------------------! 340 340 wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean 341 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice! heat sink for ocean (<0, W.m-2)341 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2) 342 342 343 343 END DO … … 381 381 CALL prt_ctl_info(' - Cell values : ') 382 382 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 383 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_itd_me : cell area :')383 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me : cell area :') 384 384 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :') 385 385 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :') … … 1093 1093 & + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 1094 1094 1095 ! in 1e-9 Joules(same as e_s)1095 ! in J/m2 (same as e_s) 1096 1096 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-fsnowrdg) & !rafting included 1097 1097 & - esrft(ji,jj)*(1.0-fsnowrft) … … 1133 1133 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux 1134 1134 1135 ! Correct dimensions to avoid big values 1136 ersw(ji,jj,jk) = ersw(ji,jj,jk) / unit_fac 1137 1138 ! Mutliply by ice volume, and divide by number of layers to get heat content in 1.e9 J 1139 ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 1140 !! MV HC 2014 1141 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) 1142 1135 ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 1143 1136 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 1144 1137 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r5055 r5064 129 129 DO jj = 1, jpj 130 130 DO ji = 1, jpi 131 rswitch 131 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 132 132 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 133 rswitch 133 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 134 134 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 135 135 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) … … 167 167 !----------------------------------------------------------------------------------------------- 168 168 !- 4.1 Compute category boundaries 169 ! Tricky trick see limitd_me.F90170 ! will be soon removed, CT171 ! hi_max(kubnd) = 99.172 169 zhbnew(:,:,:) = 0._wp 173 170 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90
r5051 r5064 114 114 CALL lbc_lnk( tmf(:,:), 'F', 1. ) ! lateral boundary conditions 115 115 116 ! !== unmasked and masked area of T-grid cell117 area(:,:) = e1t(:,:) * e2t(:,:)118 116 ! 119 117 END SUBROUTINE lim_msh -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r5053 r5064 116 116 CHARACTER (len=50) :: charout 117 117 REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta ! 118 REAL(wp) :: za, zstms , zmask! local scalars119 REAL(wp) :: zc1, zc2, zc3 120 121 REAL(wp) :: dtevp ! time step for subcycling122 REAL(wp) :: dtotel, ecc2, ecci ! square of yield ellipse eccenticity123 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars124 REAL(wp) :: zu_ice2, zv_ice1 !125 REAL(wp) :: zddc, zdtc ! delta on corners and on centre126 REAL(wp) :: zdst ! shear at the center of the grid point127 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface128 REAL(wp) :: sigma1, sigma2 ! internal ice stress118 REAL(wp) :: za, zstms ! local scalars 119 REAL(wp) :: zc1, zc2, zc3 ! ice mass 120 121 REAL(wp) :: dtevp , z1_dtevp ! time step for subcycling 122 REAL(wp) :: dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 123 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars 124 REAL(wp) :: zu_ice2, zv_ice1 ! 125 REAL(wp) :: zddc, zdtc ! delta on corners and on centre 126 REAL(wp) :: zdst ! shear at the center of the grid point 127 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface 128 REAL(wp) :: sigma1, sigma2 ! internal ice stress 129 129 130 130 REAL(wp) :: zresm ! Maximal error on ice velocity … … 137 137 REAL(wp), POINTER, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points 138 138 REAL(wp), POINTER, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays 139 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce1, v_oce1! ocean u/v component on U points140 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2 , v_oce2! ocean u/v component on V points139 REAL(wp), POINTER, DIMENSION(:,:) :: v_oce1 ! ocean u/v component on U points 140 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2 ! ocean u/v component on V points 141 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point 142 142 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses … … 156 156 157 157 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 158 CALL wrk_alloc( jpi,jpj, u_oce 1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1 )158 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 ) 159 159 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 160 160 CALL wrk_alloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) … … 228 228 ! zcorl2: Coriolis parameter on V-points 229 229 ! (ztagnx,ztagny): wind stress on U/V points 230 ! u_oce1: ocean u component on u points231 230 ! v_oce1: ocean v component on u points 232 231 ! u_oce2: ocean u component on v points 233 ! v_oce2: ocean v component on v points234 232 235 233 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==! … … 266 264 267 265 ! Mass, coriolis coeff. and currents 268 zmass1(ji,jj) = ( zt12 *zc1 + zt11*zc2 ) / (zt11+zt12+zepsi)269 zmass2(ji,jj) = ( zt22 *zc1 + zt21*zc3 ) / (zt21+zt22+zepsi)270 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) *fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) ) &266 zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 267 zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 268 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) ) & 271 269 & / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 272 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) *fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) ) &270 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) ) & 273 271 & / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 274 272 ! 275 u_oce1(ji,jj) = u_oce(ji,jj)276 v_oce2(ji,jj) = v_oce(ji,jj)277 278 273 ! Ocean has no slip boundary condition 279 v_oce1(ji,jj) = 0.5 *( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)&280 & +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj))&281 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)282 283 u_oce2(ji,jj) = 0.5 *((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)&284 & +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1))&285 & / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)274 v_oce1(ji,jj) = 0.5 * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji,jj) & 275 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 276 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 277 278 u_oce2(ji,jj) = 0.5 * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj) & 279 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 280 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 286 281 287 282 ! Wind stress at U,V-point … … 316 311 dtotel = dtevp / ( 2._wp * telast ) 317 312 #endif 313 z1_dtotel = 1._wp / ( 1._wp + dtotel ) 314 z1_dtevp = 1._wp / dtevp 318 315 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 319 316 ecc2 = ecc * ecc … … 334 331 335 332 DO jj = k_j1+1, k_jpj-1 336 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi333 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to tmi 337 334 338 335 ! … … 357 354 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 358 355 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 359 & ) / area(ji,jj)356 & ) * r1_e12t(ji,jj) 360 357 361 358 zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 362 359 & - ( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 363 & ) / area(ji,jj)360 & ) * r1_e12t(ji,jj) 364 361 365 362 ! 366 363 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 367 364 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 368 & ) / ( e1f(ji,jj) * e2f(ji,jj)) * ( 2._wp - tmf(ji,jj) ) &365 & ) * r1_e12f(ji,jj) * ( 2._wp - tmf(ji,jj) ) & 369 366 & * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1) 370 367 371 368 372 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji,jj-1) ) * e1t(ji+1,jj)&373 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) &369 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 370 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 374 371 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 375 372 376 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj) ) * e2t(ji,jj+1)&377 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) &373 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 374 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 378 375 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 379 380 END DO 381 END DO 382 CALL lbc_lnk( v_ice1, 'U', -1. ) ; CALL lbc_lnk( u_ice2, 'V', -1. ) ! lateral boundary cond. 376 END DO 377 END DO 378 CALL lbc_lnk( v_ice1 , 'U', -1. ) ; CALL lbc_lnk( u_ice2 , 'V', -1. ) ! lateral boundary cond. 383 379 384 380 DO jj = k_j1+1, k_jpj-1 … … 386 382 387 383 !- Calculate Delta at centre of grid cells 388 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) &389 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) &390 & ) / area(ji,jj)391 392 delta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst) * usecc2 )384 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 385 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) & 386 & ) * r1_e12t(ji,jj) 387 388 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 393 389 delta_i(ji,jj) = delta + creepl 394 !-Calculate stress tensor components zs1 and zs2 395 !-at centre of grid cells (see section 3.5 of CICE user's guide). 396 zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) ) & 397 & * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 398 zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) ) & 399 & / ( 1._wp + dtotel ) 400 401 END DO 402 END DO 403 404 CALL lbc_lnk( zs1(:,:), 'T', 1. ) 405 CALL lbc_lnk( zs2(:,:), 'T', 1. ) 406 407 DO jj = k_j1+1, k_jpj-1 408 DO ji = fs_2, fs_jpim1 390 409 391 !- Calculate Delta on corners 410 zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 411 & -v_ice1(ji,jj)/e1u(ji,jj) & 412 & )*e1f(ji,jj)*e1f(ji,jj) & 413 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 414 & -u_ice2(ji,jj)/e2v(ji,jj) & 415 & )*e2f(ji,jj)*e2f(ji,jj) & 416 & ) & 417 & / ( e1f(ji,jj) * e2f(ji,jj) ) 418 419 zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 420 & -v_ice1(ji,jj)/e1u(ji,jj) & 421 & )*e1f(ji,jj)*e1f(ji,jj) & 422 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 423 & -u_ice2(ji,jj)/e2v(ji,jj) & 424 & )*e2f(ji,jj)*e2f(ji,jj) & 425 & ) & 426 & / ( e1f(ji,jj) * e2f(ji,jj) ) 427 428 zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 429 430 !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 431 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * & 432 & ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) ) & 433 & / ( 1.0 + dtotel ) 434 435 END DO ! ji 436 END DO ! jj 437 438 CALL lbc_lnk( zs12(:,:), 'F', 1. ) 392 zddc = ( ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 393 & + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 394 & ) * r1_e12f(ji,jj) 395 396 zdtc = (- ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 397 & + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 398 & ) * r1_e12f(ji,jj) 399 400 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + creepl 401 402 !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 403 zs1(ji,jj) = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj) * zpresh(ji,jj) & 404 & ) * z1_dtotel 405 zs2(ji,jj) = ( zs2 (ji,jj) + dtotel * ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) & 406 & ) * z1_dtotel 407 !-Calculate stress tensor component zs12 at corners 408 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) & 409 & ) * z1_dtotel 410 411 END DO 412 END DO 413 CALL lbc_lnk( zs1 , 'T', 1. ) ; CALL lbc_lnk( zs2, 'T', 1. ) 414 CALL lbc_lnk( zs12, 'F', 1. ) 439 415 440 416 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) … … 442 418 DO ji = fs_2, fs_jpim1 443 419 !- contribution of zs1, zs2 and zs12 to zf1 444 zf1(ji,jj) = 0.5 *( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj)&445 & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj)&446 & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj)&447 & ) / ( e1u(ji,jj)*e2u(ji,jj))420 zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 421 & + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) / e2u(ji,jj) & 422 & + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) / e1u(ji,jj) & 423 & ) * r1_e12u(ji,jj) 448 424 ! contribution of zs1, zs2 and zs12 to zf2 449 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 450 & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 451 & + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 - & 452 zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 453 & ) / ( e1v(ji,jj)*e2v(ji,jj) ) 425 zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 426 & - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) / e1v(ji,jj) & 427 & + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) / e2v(ji,jj) & 428 & ) * r1_e12v(ji,jj) 454 429 END DO 455 430 END DO … … 463 438 DO jj = k_j1+1, k_jpj-1 464 439 DO ji = fs_2, fs_jpim1 465 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)466 z0 = zmass1(ji,jj) /dtevp440 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj) 441 z0 = zmass1(ji,jj) * z1_dtevp 467 442 468 443 ! SB modif because ocean has no slip boundary condition 469 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 470 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 471 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 472 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 473 (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 474 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 475 za*(u_oce1(ji,jj)) 476 zcca = z0+za 444 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji ,jj) & 445 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 446 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 447 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 448 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 449 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 450 zcca = z0 + za 477 451 zccb = zcorl1(ji,jj) 478 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+zepsi)*zmask 479 452 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 480 453 END DO 481 454 END DO … … 492 465 DO ji = fs_2, fs_jpim1 493 466 494 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)495 z0 = zmass2(ji,jj) /dtevp467 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj) 468 z0 = zmass2(ji,jj) * z1_dtevp 496 469 ! SB modif because ocean has no slip boundary condition 497 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 498 & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 499 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 500 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 501 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 502 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 503 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 504 zcca = z0+za 470 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 471 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 472 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 473 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 474 & ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 475 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 476 zcca = z0 + za 505 477 zccb = zcorl2(ji,jj) 506 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+zepsi)*zmask 507 478 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 508 479 END DO 509 480 END DO … … 520 491 DO jj = k_j1+1, k_jpj-1 521 492 DO ji = fs_2, fs_jpim1 522 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)523 z0 = zmass2(ji,jj) /dtevp493 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj) 494 z0 = zmass2(ji,jj) * z1_dtevp 524 495 ! SB modif because ocean has no slip boundary condition 525 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 526 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 527 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 528 529 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 530 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 531 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 532 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 533 zcca = z0+za 496 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 497 & +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 498 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 499 500 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 501 & ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 502 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 503 zcca = z0 + za 534 504 zccb = zcorl2(ji,jj) 535 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+zepsi)*zmask 536 505 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 537 506 END DO 538 507 END DO … … 548 517 DO jj = k_j1+1, k_jpj-1 549 518 DO ji = fs_2, fs_jpim1 550 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 551 z0 = zmass1(ji,jj)/dtevp 552 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 553 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 554 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 555 556 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 557 (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 558 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 559 za*(u_oce1(ji,jj)) 560 zcca = z0+za 519 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj) 520 z0 = zmass1(ji,jj) * z1_dtevp 521 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) & 522 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 523 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 524 525 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 526 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 527 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 528 zcca = z0 + za 561 529 zccb = zcorl1(ji,jj) 562 u_ice(ji,jj) = ( zr+zccb*zv_ice1)/(zcca+zepsi)*zmask563 END DO ! ji564 END DO ! jj530 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 531 END DO 532 END DO 565 533 566 534 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) … … 577 545 !--- Convergence test. 578 546 DO jj = k_j1+1 , k_jpj-1 579 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 580 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 581 END DO 582 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 547 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 548 END DO 549 zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) ) 583 550 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 584 551 ENDIF … … 637 604 IF ( vt_i(ji,jj) <= zvmin ) THEN 638 605 639 divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 640 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 641 & +e1v(ji,jj)*v_ice(ji,jj) & 642 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 643 & ) & 644 & / area(ji,jj) 645 646 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & 647 & -u_ice(ji-1,jj)/e2u(ji-1,jj) & 648 & )*e2t(ji,jj)*e2t(ji,jj) & 649 & -( v_ice(ji,jj)/e1v(ji,jj) & 650 & -v_ice(ji,jj-1)/e1v(ji,jj-1) & 651 & )*e1t(ji,jj)*e1t(ji,jj) & 652 & ) & 653 & / area(ji,jj) 606 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj ) * u_ice(ji-1,jj ) & 607 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji ,jj-1) * v_ice(ji ,jj-1) & 608 & ) * r1_e12t(ji,jj) 609 610 zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 611 & -( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 612 & ) * r1_e12t(ji,jj) 654 613 ! 655 614 ! SB modif because ocean has no slip boundary condition 656 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) & 657 & - u_ice(ji,jj) / e1u(ji,jj) ) & 658 & * e1f(ji,jj) * e1f(ji,jj) & 659 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) & 660 & - v_ice(ji,jj) / e2v(ji,jj) ) & 661 & * e2f(ji,jj) * e2f(ji,jj) ) & 662 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 663 & * tmi(ji,jj) * tmi(ji,jj+1) & 664 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 665 666 zdst = ( e2u( ji , jj ) * v_ice1(ji ,jj ) & 667 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj ) & 668 & + e1v( ji , jj ) * u_ice2(ji ,jj ) & 669 & - e1v( ji , jj-1 ) * u_ice2(ji ,jj-1) ) / area(ji,jj) 670 671 delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 615 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 616 & +( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 617 & ) * r1_e12f(ji,jj) * ( 2.0 - tmf(ji,jj) ) & 618 & * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1) 619 620 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 621 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) ) * r1_e12t(ji,jj) 622 623 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 672 624 delta_i(ji,jj) = delta + creepl 673 625 674 626 ENDIF 675 676 END DO !jj 677 END DO !ji 627 END DO 628 END DO 678 629 ! 679 630 !------------------------------------------------------------------------------! … … 684 635 DO jj = k_j1+1, k_jpj-1 685 636 DO ji = fs_2, fs_jpim1 686 zdst= ( e2u( ji , jj ) * v_ice1(ji,jj) & 687 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 688 & + e1v( ji , jj ) * u_ice2(ji,jj) & 689 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 637 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 638 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 690 639 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 691 640 END DO … … 741 690 ! 742 691 CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 743 CALL wrk_dealloc( jpi,jpj, u_oce 1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1 )692 CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 ) 744 693 CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 745 694 CALL wrk_dealloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r5055 r5064 26 26 USE phycst ! physical constants 27 27 USE dom_oce ! ocean domain 28 USE dom_ice, ONLY : tms , area28 USE dom_ice, ONLY : tms 29 29 USE ice ! LIM sea-ice variables 30 30 USE sbc_ice ! Surface boundary condition: sea-ice fields … … 99 99 !! Tartinville et al. 2001 Ocean Modelling, 3, 95-108. 100 100 !! These refs are now obsolete since everything has been revised 101 !! The ref should be Rousset et al., 2015 ?101 !! The ref should be Rousset et al., 2015 102 102 !!--------------------------------------------------------------------- 103 103 INTEGER, INTENT(in) :: kt ! number of iteration 104 !105 104 INTEGER :: ji, jj, jl, jk ! dummy loop indices 106 ! 107 REAL(wp) :: zemp ! local scalars 105 REAL(wp) :: zemp ! local scalars 108 106 REAL(wp) :: zf_mass ! Heat flux associated with mass exchange ice->ocean (W.m-2) 109 107 REAL(wp) :: zfcm1 ! New solar flux received by the ocean … … 321 319 !! ** input : Namelist namicedia 322 320 !!------------------------------------------------------------------- 323 REAL(wp) :: zsum, zarea324 !325 321 INTEGER :: ji, jj, jk ! dummy loop indices 326 322 REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5055 r5064 109 109 ! 1.2) Heat content 110 110 !-------------------- 111 ! Change the units of heat content; from global units to J.m3111 ! Change the units of heat content; from J/m2 to J/m3 112 112 DO jl = 1, jpl 113 113 DO jk = 1, nlay_i … … 115 115 DO ji = 1, jpi 116 116 !0 if no ice and 1 if yes 117 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi 10 ) )117 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 ) ) 118 118 !Energy of melting q(S,T) [J.m-3] 119 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 ) 120 !convert units ! very important that this line is here 121 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 119 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i ) 122 120 END DO 123 121 END DO … … 127 125 DO ji = 1, jpi 128 126 !0 if no ice and 1 if yes 129 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi 10 ) )127 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi20 ) ) 130 128 !Energy of melting q(S,T) [J.m-3] 131 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 ) 132 !convert units ! very important that this line is here 133 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac 129 e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s ) 134 130 END DO 135 131 END DO … … 453 449 ! Ice heat content 454 450 !------------------------ 455 ! Enthalpies are global variables we have to readjust the units (heat content in J oules)451 ! Enthalpies are global variables we have to readjust the units (heat content in J/m2) 456 452 DO jl = 1, jpl 457 453 DO jk = 1, nlay_i 458 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a rea(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) / ( unit_fac * REAL( nlay_i ))454 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) / REAL( nlay_i ) 459 455 END DO 460 456 END DO … … 463 459 ! Snow heat content 464 460 !------------------------ 465 ! Enthalpies are global variables we have to readjust the units (heat content in J oules)461 ! Enthalpies are global variables we have to readjust the units (heat content in J/m2) 466 462 DO jl = 1, jpl 467 463 DO jk = 1, nlay_s 468 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a rea(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) / ( unit_fac * REAL( nlay_s ))464 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) / REAL( nlay_s ) 469 465 END DO 470 466 END DO … … 489 485 CALL prt_ctl_info(' - Cell values : ') 490 486 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 491 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_thd : cell area :')487 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_thd : cell area :') 492 488 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd : at_i :') 493 489 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd : vt_i :') … … 520 516 CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 521 517 518 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 522 519 !------------------------------------------------------------------------------| 523 520 ! 1) Transport of ice between thickness categories. | 524 521 !------------------------------------------------------------------------------| 525 ! Given thermodynamic growth rates, transport ice between 526 ! thickness categories. 522 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 523 524 ! Given thermodynamic growth rates, transport ice between thickness categories. 527 525 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 528 526 ! … … 530 528 CALL lim_var_agg(1) 531 529 530 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 531 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 532 532 !------------------------------------------------------------------------------| 533 533 ! 3) Add frazil ice growing in leads. … … 536 536 CALL lim_var_glo2eqv ! only for info 537 537 538 ! conservation test 539 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 540 538 541 IF(ln_ctl) THEN ! Control print 539 542 CALL prt_ctl_info(' ') 540 543 CALL prt_ctl_info(' - Cell values : ') 541 544 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 542 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_itd_th : cell area :')545 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_th : cell area :') 543 546 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th : at_i :') 544 547 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th : vt_i :') … … 568 571 ENDIF 569 572 ! 570 ! conservation test571 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)572 573 573 IF( nn_timing == 1 ) CALL timing_stop('limthd') 574 574 … … 614 614 !! ( dA = A/2h dh ) 615 615 !!----------------------------------------------------------------------- 616 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 617 INTEGER :: ji ! dummy loop indices 618 REAL(wp) :: zhi_bef ! ice thickness before thermo 619 REAL(wp) :: zdh_mel ! net melting 616 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 617 INTEGER :: ji ! dummy loop indices 618 REAL(wp) :: zhi_bef ! ice thickness before thermo 619 REAL(wp) :: zdh_mel, zda_mel ! net melting 620 REAL(wp) :: zv ! ice volume 620 621 621 622 DO ji = kideb, kiut 622 zhi_bef = ht_i_1d(ji) - ( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 623 zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 624 IF( zdh_mel < 0._wp ) & 625 & a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) ) 623 zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 624 IF( zdh_mel < 0._wp ) THEN 625 zv = a_i_1d(ji) * ht_i_1d(ji) 626 ! lateral melting = concentration change 627 zhi_bef = ht_i_1d(ji) - zdh_mel 628 zda_mel = a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) 629 a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + zda_mel ) 630 ! adjust thickness 631 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - a_i_1d(ji) + epsi20 ) ) 632 ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 633 ! adjust thickness 634 at_i_1d(ji) = a_i_1d(ji) 635 END IF 626 636 END DO 627 at_i_1d(:) = a_i_1d(:)628 637 629 638 END SUBROUTINE lim_thd_lam … … 665 674 IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 666 675 nn_monocat = 0 667 WRITE(numout, *) 'nn_monocat must be 0 in multi-category case '676 IF(lwp) WRITE(numout, *) ' nn_monocat must be 0 in multi-category case ' 668 677 ENDIF 669 678 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5055 r5064 157 157 REAL(wp) :: zh_thres ! thickness thres. for G(h) computation 158 158 REAL(wp) :: zhe ! dummy factor 159 REAL(wp) :: zswitch ! dummy switch160 159 REAL(wp) :: zkimean ! mean sea ice thermal conductivity 161 160 REAL(wp) :: zfac ! dummy factor … … 363 362 ! Fichefet and Morales Maqueda, JGR 1997 364 363 365 zghe(:) = 0._wp364 zghe(:) = 1._wp 366 365 367 366 SELECT CASE ( nn_monocat ) 368 369 CASE (0,2,4)370 371 zghe(kideb:kiut) = 1._wp372 367 373 368 CASE (1,3) ! LIM3 … … 388 383 389 384 ! G(he) 390 zswitch = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) ) ! =0 if zhe < zh_thres, if >391 zghe(ji) = ( 1. 0 - zswitch ) + zswitch * ( 0.5 + 0.5 *LOG( 2.*zhe / zepsilon ) )385 rswitch = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) ) ! =0 if zhe < zh_thres, if > 386 zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2.*zhe / zepsilon ) ) 392 387 393 388 ! Impose G(he) < 2. 394 zghe(ji) = MIN( zghe(ji), 2. 0)389 zghe(ji) = MIN( zghe(ji), 2._wp ) 395 390 396 391 END DO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r5055 r5064 132 132 DO jk1 = 1, nlay_i 133 133 DO ji = kideb, kiut 134 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi 10 ) )135 qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi 10 )134 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi20 ) ) 135 qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 136 136 ENDDO 137 137 ENDDO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r5055 r5064 130 130 DO ji = 1, jpi 131 131 !Energy of melting q(S,T) [J.m-3] 132 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 133 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) & 134 & / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i, wp ) 135 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 132 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi20 ) ) !0 if no ice 133 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp ) 136 134 END DO 137 135 END DO … … 526 524 DO jj = 1, jpj 527 525 DO ji = 1, jpi 528 ! heat content in J oules529 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac)526 ! heat content in J/m2 527 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / REAL( nlay_i ,wp ) 530 528 END DO 531 529 END DO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r5059 r5064 95 95 ENDIF 96 96 97 zsm(:,:) = area(:,:)97 zsm(:,:) = e12t(:,:) 98 98 99 99 ! !-------------------------------------! … … 153 153 ! transported fields 154 154 !------------------------- 155 zs0ow(:,:,1) = ato_i(:,:) * area(:,:) ! Open water area 156 DO jl = 1, jpl 157 zs0sn (:,:,jl) = v_s (:,:,jl) * area(:,:) ! Snow volume 158 zs0ice(:,:,jl) = v_i (:,:,jl) * area(:,:) ! Ice volume 159 zs0a (:,:,jl) = a_i (:,:,jl) * area(:,:) ! Ice area 160 zs0sm (:,:,jl) = smv_i(:,:,jl) * area(:,:) ! Salt content 161 zs0oi (:,:,jl) = oa_i (:,:,jl) * area(:,:) ! Age content 162 zs0c0 (:,:,jl) = e_s (:,:,1,jl) ! Snow heat content 163 zs0e (:,:,:,jl) = e_i (:,:,:,jl) ! Ice heat content 155 zs0ow(:,:,1) = ato_i(:,:) * e12t(:,:) ! Open water area 156 DO jl = 1, jpl 157 zs0sn (:,:,jl) = v_s (:,:,jl) * e12t(:,:) ! Snow volume 158 zs0ice(:,:,jl) = v_i (:,:,jl) * e12t(:,:) ! Ice volume 159 zs0a (:,:,jl) = a_i (:,:,jl) * e12t(:,:) ! Ice area 160 zs0sm (:,:,jl) = smv_i(:,:,jl) * e12t(:,:) ! Salt content 161 zs0oi (:,:,jl) = oa_i (:,:,jl) * e12t(:,:) ! Age content 162 zs0c0 (:,:,jl) = e_s (:,:,1,jl) * e12t(:,:) ! Snow heat content 163 DO jk = 1, nlay_i 164 zs0e (:,:,jk,jl) = e_i (:,:,jk,jl) * e12t(:,:) ! Ice heat content 165 END DO 164 166 END DO 165 167 … … 253 255 ! Recover the properties from their contents 254 256 !------------------------------------------- 255 ato_i(:,:) = zs0ow(:,:,1) / area(:,:) 256 DO jl = 1, jpl 257 v_i (:,:,jl) = zs0ice(:,:,jl) / area(:,:) 258 v_s (:,:,jl) = zs0sn (:,:,jl) / area(:,:) 259 smv_i(:,:,jl) = zs0sm (:,:,jl) / area(:,:) 260 oa_i (:,:,jl) = zs0oi (:,:,jl) / area(:,:) 261 a_i (:,:,jl) = zs0a (:,:,jl) / area(:,:) 262 e_s (:,:,1,jl) = zs0c0 (:,:,jl) 263 e_i (:,:,:,jl) = zs0e (:,:,:,jl) 257 ato_i(:,:) = zs0ow(:,:,1) * r1_e12t(:,:) 258 DO jl = 1, jpl 259 v_i (:,:,jl) = zs0ice(:,:,jl) * r1_e12t(:,:) 260 v_s (:,:,jl) = zs0sn (:,:,jl) * r1_e12t(:,:) 261 smv_i(:,:,jl) = zs0sm (:,:,jl) * r1_e12t(:,:) 262 oa_i (:,:,jl) = zs0oi (:,:,jl) * r1_e12t(:,:) 263 a_i (:,:,jl) = zs0a (:,:,jl) * r1_e12t(:,:) 264 e_s (:,:,1,jl) = zs0c0 (:,:,jl) * r1_e12t(:,:) 265 DO jk = 1, nlay_i 266 e_i (:,:,jk,jl) = zs0e (:,:,jk,jl) * r1_e12t(:,:) 267 END DO 264 268 END DO 265 269 … … 380 384 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 381 385 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 382 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 <0383 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 <0386 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 387 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0 384 388 ENDIF 385 389 … … 404 408 DO jj = 1, jpj 405 409 DO ji = 1, jpi 406 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac* r1_rdtice407 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac* r1_rdtice410 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 411 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 408 412 409 413 diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r5059 r5064 16 16 USE sbc_ice ! Surface boundary condition: ice fields 17 17 USE dom_ice 18 USE dom_oce 18 19 USE phycst ! physical constants 19 20 USE ice … … 148 149 diag_heat_dhc(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) + & 149 150 & SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) & 150 & ) * unit_fac * r1_rdtice / area(ji,jj)151 & ) * r1_rdtice 151 152 END DO 152 153 END DO … … 159 160 CALL prt_ctl_info(' - Cell values : ') 160 161 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 161 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_update1 : cell area :')162 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_update1 : cell area :') 162 163 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_update1 : at_i :') 163 164 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_update1 : vt_i :') … … 183 184 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_update1 : v_s : ') 184 185 CALL prt_ctl(tab2d_1=v_s_b (:,:,jl) , clinfo1= ' lim_update1 : v_s_b : ') 185 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl) /1.0e9, clinfo1= ' lim_update1 : e_i1 : ')186 CALL prt_ctl(tab2d_1=e_i_b (:,:,1,jl) /1.0e9, clinfo1= ' lim_update1 : e_i1_b : ')187 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl) /1.0e9, clinfo1= ' lim_update1 : e_i2 : ')188 CALL prt_ctl(tab2d_1=e_i_b (:,:,2,jl) /1.0e9, clinfo1= ' lim_update1 : e_i2_b : ')186 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl) , clinfo1= ' lim_update1 : e_i1 : ') 187 CALL prt_ctl(tab2d_1=e_i_b (:,:,1,jl) , clinfo1= ' lim_update1 : e_i1_b : ') 188 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl) , clinfo1= ' lim_update1 : e_i2 : ') 189 CALL prt_ctl(tab2d_1=e_i_b (:,:,2,jl) , clinfo1= ' lim_update1 : e_i2_b : ') 189 190 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_update1 : e_snow : ') 190 191 CALL prt_ctl(tab2d_1=e_s_b (:,:,1,jl) , clinfo1= ' lim_update1 : e_snow_b : ') -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r5059 r5064 16 16 USE sbc_ice ! Surface boundary condition: ice fields 17 17 USE dom_ice 18 USE dom_oce 18 19 USE phycst ! physical constants 19 20 USE ice … … 190 191 & ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) + & 191 192 & SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) & 192 & ) * unit_fac * r1_rdtice / area(ji,jj)193 & ) * r1_rdtice 193 194 END DO 194 195 END DO … … 203 204 CALL prt_ctl_info(' - Cell values : ') 204 205 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 205 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_update2 : cell area :')206 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_update2 : cell area :') 206 207 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_update2 : at_i :') 207 208 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_update2 : vt_i :') … … 227 228 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_update2 : v_s : ') 228 229 CALL prt_ctl(tab2d_1=v_s_b (:,:,jl) , clinfo1= ' lim_update2 : v_s_b : ') 229 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl) /1.0e9, clinfo1= ' lim_update2 : e_i1 : ')230 CALL prt_ctl(tab2d_1=e_i_b (:,:,1,jl) /1.0e9, clinfo1= ' lim_update2 : e_i1_b : ')231 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl) /1.0e9, clinfo1= ' lim_update2 : e_i2 : ')232 CALL prt_ctl(tab2d_1=e_i_b (:,:,2,jl) /1.0e9, clinfo1= ' lim_update2 : e_i2_b : ')230 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl) , clinfo1= ' lim_update2 : e_i1 : ') 231 CALL prt_ctl(tab2d_1=e_i_b (:,:,1,jl) , clinfo1= ' lim_update2 : e_i1_b : ') 232 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl) , clinfo1= ' lim_update2 : e_i2 : ') 233 CALL prt_ctl(tab2d_1=e_i_b (:,:,2,jl) , clinfo1= ' lim_update2 : e_i2_b : ') 233 234 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_update2 : e_snow : ') 234 235 CALL prt_ctl(tab2d_1=e_s_b (:,:,1,jl) , clinfo1= ' lim_update2 : e_snow_b : ') -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r5055 r5064 190 190 DO ji = 1, jpi 191 191 ! ! Energy of melting q(S,T) [J.m-3] 192 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) ! rswitch = 0 if no ice and 1 if yes 193 zq_i = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 194 zq_i = zq_i * unit_fac ! convert units 192 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 193 zq_i = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 195 194 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt ! Ice layer melt temperature 196 195 ! … … 216 215 DO ji = 1, jpi 217 216 !Energy of melting q(S,T) [J.m-3] 218 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) ) ! rswitch = 0 if no ice and 1 if yes 219 zq_s = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp) 220 zq_s = zq_s * unit_fac ! convert units 217 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 218 zq_s = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp) 221 219 ! 222 220 t_s(ji,jj,jk,jl) = rtt + rswitch * ( - zfac1 * zq_s + zfac2 ) … … 541 539 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rtt * ( 1._wp - rswitch ) 542 540 ! update exchanges with ocean 543 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 <0541 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 544 542 END DO 545 543 END DO … … 579 577 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 580 578 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 581 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 <0579 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 582 580 END DO 583 581 END DO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90
r5055 r5064 229 229 ! Snow energy of melting 230 230 e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 231 ! Change dimensions 232 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 233 ! Multiply by volume, so that heat content in 10^9 Joules 234 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 231 ! Multiply by volume, so that heat content in J/m2 232 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) / nlay_s 235 233 END DO 236 234 DO jk = 1, nlay_i … … 241 239 + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 242 240 - rcp * ( ztmelts - rtt ) ) 243 ! Correct dimensions to avoid big values 244 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 245 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 246 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i 241 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 242 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i 247 243 END DO 248 244
Note: See TracChangeset
for help on using the changeset viewer.
