Changeset 4099 for branches/2013/dev_r4028_CNRS_LIM3
- Timestamp:
- 2013-10-22T14:07:21+02:00 (10 years ago)
- Location:
- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM
- Files:
-
- 2 added
- 2 deleted
- 17 edited
-
CONFIG/ORCA2_LIM3/MY_SRC (added)
-
CONFIG/ORCA2_LIM3/cpp_ORCA2_LIM3.fcm (added)
-
CONFIG/cfg.txt (modified) (1 diff)
-
NEMO/LIM_SRC_3/ice.F90 (modified) (2 diffs)
-
NEMO/LIM_SRC_3/iceini.F90 (modified) (2 diffs)
-
NEMO/LIM_SRC_3/limcat_1D.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limcat_2D.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limdia.F90 (deleted)
-
NEMO/LIM_SRC_3/limistate.F90 (modified) (10 diffs)
-
NEMO/LIM_SRC_3/limrhg.F90 (modified) (2 diffs)
-
NEMO/LIM_SRC_3/limrst.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limsbc.F90 (modified) (3 diffs)
-
NEMO/LIM_SRC_3/limthd_dh.F90 (modified) (3 diffs)
-
NEMO/LIM_SRC_3/limupdate.F90 (deleted)
-
NEMO/LIM_SRC_3/limupdate1.F90 (modified) (3 diffs)
-
NEMO/LIM_SRC_3/limupdate2.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/SBC/sbc_ice.F90 (modified) (1 diff)
-
NEMO/OPA_SRC/SBC/sbccpl.F90 (modified) (3 diffs)
-
NEMO/OPA_SRC/SBC/sbcice_if.F90 (modified) (1 diff)
-
NEMO/OPA_SRC/SBC/sbcice_lim.F90 (modified) (8 diffs)
-
NEMO/OPA_SRC/SBC/sbcmod.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/cfg.txt
r3905 r4099 1 AMM12 OPA_SRC2 1 GYRE OPA_SRC 3 2 GYRE_BFM OPA_SRC TOP_SRC 4 3 GYRE_PISCES OPA_SRC TOP_SRC 5 ORCA2_LIM3 OPA_SRC LIM_SRC_36 4 ORCA2_SAS_LIM OPA_SRC SAS_SRC LIM_SRC_2 NST_SRC 7 5 ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 8 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC9 6 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 10 7 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 8 AMM12 OPA_SRC 9 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 10 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r4045 r4099 395 395 CHARACTER(len=32) , PUBLIC :: cn_icerst_out = "restart_ice" !: suffix of ice restart name (output) 396 396 LOGICAL , PUBLIC :: ln_limdyn = .TRUE. !: flag for ice dynamics (T) or not (F) 397 LOGICAL , PUBLIC :: ln_nicep = . TRUE. !: flag for sea-ice points output (T) or not (F)397 LOGICAL , PUBLIC :: ln_nicep = .FALSE. !: flag for sea-ice points output (T) or not (F) 398 398 REAL(wp) , PUBLIC :: cai = 1.40e-3 !: atmospheric drag over sea ice 399 399 REAL(wp) , PUBLIC :: cao = 1.00e-3 !: atmospheric drag over ocean … … 404 404 !!-------------------------------------------------------------------------- 405 405 !! Check if everything down here is necessary 406 LOGICAL , PUBLIC :: ln_limdiahsb = .TRUE. !: flag for ice diag (T) or not (F) 406 LOGICAL , PUBLIC :: ln_limdiahsb = .FALSE. !: flag for ice diag (T) or not (F) 407 LOGICAL , PUBLIC :: ln_limdiaout = .FALSE. !: flag for ice diag (T) or not (F) 407 408 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_newice !: volume of ice formed in the leads 408 409 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dv_dt_thd !: thermodynamic growth rates -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
r4045 r4099 126 126 !! ** input : Namelist namicerun 127 127 !!------------------------------------------------------------------- 128 NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb 128 NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb, ln_limdiaout 129 129 !!------------------------------------------------------------------- 130 130 ! … … 147 147 WRITE(numout,*) ' Several ice points in the ice or not in ocean.output = ', ln_nicep 148 148 WRITE(numout,*) ' Diagnose heat/salt budget or not ln_limdiahsb = ', ln_limdiahsb 149 WRITE(numout,*) ' Output heat/salt budget or not ln_limdiaout = ', ln_limdiaout 149 150 ENDIF 150 151 ! -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90
r4072 r4099 154 154 !============================ 155 155 ! Check if tests have passed (i.e. volume conservation...) 156 IF ( ztests /= 4 ) THEN157 WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '158 WRITE(numout,*) ' !! ALERT categories distribution !!'159 WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '160 WRITE(numout,*) ' *** ztests is not equal to 4 '161 WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4162 WRITE(numout,*) 'i_fill=',i_fill163 WRITE(numout,*) 'zai(ji)=',zai(ji)164 WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:)165 ENDIF156 !IF ( ztests /= 4 ) THEN 157 ! WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 158 ! WRITE(numout,*) ' !! ALERT categories distribution !!' 159 ! WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 160 ! WRITE(numout,*) ' *** ztests is not equal to 4 ' 161 ! WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4 162 ! WRITE(numout,*) 'i_fill=',i_fill 163 ! WRITE(numout,*) 'zai(ji)=',zai(ji) 164 ! WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:) 165 !ENDIF 166 166 167 167 END DO ! i loop -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_2D.F90
r4072 r4099 155 155 !============================ 156 156 ! Check if tests have passed (i.e. volume conservation...) 157 IF ( ztests .NE. 4 ) THEN158 WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '159 WRITE(numout,*) ' !! ALERT categories distribution !!'160 WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '161 WRITE(numout,*) ' *** ztests is not equal to 4 '162 WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4163 ENDIF157 !IF ( ztests .NE. 4 ) THEN 158 ! WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 159 ! WRITE(numout,*) ' !! ALERT categories distribution !!' 160 ! WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 161 ! WRITE(numout,*) ' *** ztests is not equal to 4 ' 162 ! WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4 163 !ENDIF 164 164 165 165 END DO ! i loop -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r4045 r4099 19 19 USE dom_oce ! ocean domain 20 20 USE sbc_oce ! Surface boundary condition: ocean fields 21 USE sbc_ice ! Surface boundary condition: ice fields 21 22 USE eosbn2 ! equation of state 22 23 USE ice ! sea-ice variables … … 27 28 USE lbclnk ! lateral boundary condition - MPP exchanges 28 29 USE lib_mpp ! MPP library 30 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 29 31 USE wrk_nemo ! work arrays 30 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)31 32 32 33 IMPLICIT NONE … … 35 36 PUBLIC lim_istate ! routine called by lim_init.F90 36 37 37 ! !!** init namelist (namiceini) ** 38 REAL(wp) :: ttest = 2.0_wp ! threshold water temperature for initial sea ice 39 REAL(wp) :: hninn = 0.5_wp ! initial snow thickness in the north 40 REAL(wp) :: hginn_u = 2.5_wp ! initial ice thickness in the north 41 REAL(wp) :: aginn_u = 0.7_wp ! initial leads area in the north 42 REAL(wp) :: hginn_d = 5.0_wp ! initial ice thickness in the north 43 REAL(wp) :: aginn_d = 0.25_wp ! initial leads area in the north 44 REAL(wp) :: hnins = 0.1_wp ! initial snow thickness in the south 45 REAL(wp) :: hgins_u = 1.0_wp ! initial ice thickness in the south 46 REAL(wp) :: agins_u = 0.7_wp ! initial leads area in the south 47 REAL(wp) :: hgins_d = 2.0_wp ! initial ice thickness in the south 48 REAL(wp) :: agins_d = 0.2_wp ! initial leads area in the south 49 REAL(wp) :: sinn = 6.301_wp ! initial salinity 50 REAL(wp) :: sins = 6.301_wp ! 51 52 !!---------------------------------------------------------------------- 53 !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 38 !! * Module variables 39 REAL(wp) :: & !!! ** init namelist (namiceini) ** 40 ttest = 2.0 , & ! threshold water temperature for initial sea ice 41 hninn = 0.5 , & ! initial snow thickness in the north 42 hnins = 0.1 , & ! initial snow thickness in the south 43 hginn = 2.5 , & ! initial ice thickness in the north 44 hgins = 1.0 , & ! initial ice thickness in the south 45 aginn = 0.7 , & ! initial leads area in the north 46 agins = 0.7 , & ! initial leads area in the south 47 sinn = 6.301 , & ! initial salinity 48 sins = 6.301 49 50 !!---------------------------------------------------------------------- 51 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) 54 52 !! $Id$ 55 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 56 !!---------------------------------------------------------------------- 53 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 54 !!---------------------------------------------------------------------- 55 57 56 CONTAINS 58 57 … … 63 62 !! ** Purpose : defined the sea-ice initial state 64 63 !! 65 !! ** Method : restart from a state defined in a binary file 66 !! or from arbitrary sea-ice conditions 67 !!------------------------------------------------------------------- 68 INTEGER :: ji, jj, jk, jl ! dummy loop indices 69 REAL(wp) :: zeps6, zeps, ztmelts, epsi06 ! local scalars 70 REAL(wp) :: zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs 71 REAL(wp), POINTER, DIMENSION(:) :: zgfactorn, zhin 72 REAL(wp), POINTER, DIMENSION(:) :: zgfactors, zhis 73 REAL(wp), POINTER, DIMENSION(:,:) :: zidto ! ice indicator 74 !-------------------------------------------------------------------- 75 76 CALL wrk_alloc( jpm, zgfactorn, zgfactors, zhin, zhis ) 64 !! ** Method : 65 !! This routine will put some ice where ocean 66 !! is at the freezing point, then fill in ice 67 !! state variables using prescribed initial 68 !! values in the namelist 69 !! 70 !! ** Steps : 71 !! 1) Read namelist 72 !! 2) Basal temperature; ice and hemisphere masks 73 !! 3) Fill in the ice thickness distribution using gaussian 74 !! 4) Fill in space-dependent arrays for state variables 75 !! 5) Diagnostic arrays 76 !! 6) Lateral boundary conditions 77 !! 78 !! History : 79 !! 2.0 ! 01-04 (C. Ethe, G. Madec) Original code 80 !! 3.0 ! 2007 (M. Vancoppenolle) Rewrite for ice cats 81 !! 4.0 ! 09-11 (M. Vancoppenolle) Enhanced version for ice cats 82 !!-------------------------------------------------------------------- 83 84 !! * Local variables 85 INTEGER :: ji, jj, jk, jl ! dummy loop indices 86 REAL(wp) :: epsi06, epsi20, ztmelts 87 INTEGER :: i_hemis, i_fill, jl0 88 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 89 REAL(wp), POINTER, DIMENSION(:) :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini 90 REAL(wp), POINTER, DIMENSION(:,:) :: zht_i_ini, za_i_ini, zv_i_ini 91 REAL(wp), POINTER, DIMENSION(:,:) :: zidto ! ice indicator 92 INTEGER, POINTER, DIMENSION(:,:) :: zhemis ! hemispheric index 93 !-------------------------------------------------------------------- 94 77 95 CALL wrk_alloc( jpi, jpj, zidto ) 78 79 !-------------------------------------------------------------------- 80 ! 1) Preliminary things 81 !-------------------------------------------------------------------- 82 epsi06 = 1.e-6_wp 96 CALL wrk_alloc( jpi, jpj, zhemis ) 97 CALL wrk_alloc( jpl, 2, zht_i_ini, za_i_ini, zv_i_ini ) 98 CALL wrk_alloc( 2, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 99 100 epsi06 = 1.0e-6 101 epsi20 = 1.0e-20 102 IF(lwp) WRITE(numout,*) 103 IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' 104 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 105 106 !-------------------------------------------------------------------- 107 ! 1) Read namelist 108 !-------------------------------------------------------------------- 83 109 84 110 CALL lim_istate_init ! reading the initials parameters of the ice … … 89 115 90 116 !-------------------------------------------------------------------- 91 ! 2) Ice initialization (hi,hs,frld,t_su,sm_i,t_i,t_s) | 92 !-------------------------------------------------------------------- 93 94 IF(lwp) WRITE(numout,*) 95 IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' 96 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 97 117 ! 2) Basal temperature, ice mask and hemispheric index 118 !-------------------------------------------------------------------- 119 120 ! Basal temperature is set to the freezing point of seawater in Celsius 98 121 t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius] 99 122 100 123 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 101 124 DO ji = 1, jpi 102 IF( tsn(ji,jj,1,jp_tem) - t_bo(ji,jj) >= ttest ) THEN ; zidto(ji,jj) = 0. e0! no ice103 ELSE ; zidto(ji,jj) = 1. e0! ice125 IF( tsn(ji,jj,1,jp_tem) - t_bo(ji,jj) >= ttest ) THEN ; zidto(ji,jj) = 0._wp ! no ice 126 ELSE ; zidto(ji,jj) = 1._wp ! ice 104 127 ENDIF 105 128 END DO 106 129 END DO 107 130 108 t_bo(:,:) = t_bo(:,:) + rt0 ! t_bo converted from Celsius to Kelvin (rt0 over land) 109 110 ! constants for heat contents 111 zeps = 1.e-20_wp 112 zeps6 = 1.e-06_wp 113 114 ! zgfactor for initial ice distribution 115 zgfactorn(:) = 0._wp 116 zgfactors(:) = 0._wp 117 118 ! first ice type 119 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 120 zhin (1) = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 121 zgfactorn(1) = zgfactorn(1) + exp(-(zhin(1)-hginn_u)*(zhin(1)-hginn_u) * 0.5_wp ) 122 zhis (1) = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 123 zgfactors(1) = zgfactors(1) + exp(-(zhis(1)-hgins_u)*(zhis(1)-hgins_u) * 0.5_wp ) 124 END DO ! jl 125 zgfactorn(1) = aginn_u / zgfactorn(1) 126 zgfactors(1) = agins_u / zgfactors(1) 127 128 ! ------------- 129 ! new distribution, polynom of second order, conserving area and volume 130 zh1 = 0._wp 131 zh2 = 0._wp 132 zh3 = 0._wp 133 DO jl = 1, jpl 134 zh = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 135 zh1 = zh1 + zh 136 zh2 = zh2 + zh * zh 137 zh3 = zh3 + zh * zh * zh 138 END DO 139 IF(lwp) WRITE(numout,*) ' zh1 : ', zh1 140 IF(lwp) WRITE(numout,*) ' zh2 : ', zh2 141 IF(lwp) WRITE(numout,*) ' zh3 : ', zh3 142 143 zvol = aginn_u * hginn_u 144 zare = aginn_u 145 IF( jpl >= 2 ) THEN 146 zbn = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3) 147 zan = ( zare - zbn*zh1 ) / zh2 148 ENDIF 149 150 IF(lwp) WRITE(numout,*) ' zvol: ', zvol 151 IF(lwp) WRITE(numout,*) ' zare: ', zare 152 IF(lwp) WRITE(numout,*) ' zbn : ', zbn 153 IF(lwp) WRITE(numout,*) ' zan : ', zan 154 155 zvol = agins_u * hgins_u 156 zare = agins_u 157 IF( jpl >= 2 ) THEN 158 zbs = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3) 159 zas = ( zare - zbs*zh1 ) / zh2 160 ENDIF 161 162 IF(lwp) WRITE(numout,*) ' zvol: ', zvol 163 IF(lwp) WRITE(numout,*) ' zare: ', zare 164 IF(lwp) WRITE(numout,*) ' zbn : ', zbn 165 IF(lwp) WRITE(numout,*) ' zan : ', zan 166 167 !end of new lines 168 ! ------------- 169 !!! 170 ! retour a LIMA_MEC 171 ! ! second ice type 172 ! zdummy = hi_max(ice_cat_bounds(2,1)-1) 173 ! hi_max(ice_cat_bounds(2,1)-1) = 0.0 174 175 ! ! here to change !!!! 176 ! jm = 2 177 ! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 178 ! zhin (2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 179 ! zhin (2) = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm ) + & 180 ! hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm) ) / 2.0 181 ! zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0) 182 ! zhis (2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 183 ! zhis (2) = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm ) + & 184 ! hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm) ) / 2.0 185 ! zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0) 186 ! END DO ! jl 187 ! zgfactorn(2) = aginn_d / zgfactorn(2) 188 ! zgfactors(2) = agins_d / zgfactors(2) 189 ! hi_max(ice_cat_bounds(2,1)-1) = zdummy 190 ! END retour a LIMA_MEC 191 !!! 192 193 !!gm optimisation : loop over the ice categories inside the ji, jj loop !!! 194 131 t_bo(:,:) = t_bo(:,:) + rt0 ! conversion to Kelvin 132 133 ! Hemispheric index 134 ! MV 2011 new initialization 195 135 DO jj = 1, jpj 196 136 DO ji = 1, jpi 197 198 !--- Northern hemisphere199 !----------------------------------------------------------------200 137 IF( fcor(ji,jj) >= 0._wp ) THEN 201 202 !----------------------- 203 ! Ice area / thickness 204 !----------------------- 205 206 IF ( jpl .EQ. 1) THEN ! one category 207 208 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 209 a_i(ji,jj,jl) = zidto(ji,jj) * aginn_u 210 ht_i(ji,jj,jl) = zidto(ji,jj) * hginn_u 211 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 212 END DO 213 214 ELSE ! several categories 215 216 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 217 zhin(1) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 218 a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* & 219 (zhin(1)-hginn_u)/2.0) , epsi06) 220 ! new line 221 a_i(ji,jj,jl) = zidto(ji,jj) * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) ) 222 ht_i(ji,jj,jl) = zidto(ji,jj) * zhin(1) 223 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 224 END DO 225 226 ENDIF 227 228 229 !!! 230 ! retour a LIMA_MEC 231 ! !ridged ice 232 ! zdummy = hi_max(ice_cat_bounds(2,1)-1) 233 ! hi_max(ice_cat_bounds(2,1)-1) = 0.0 234 ! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories 235 ! zhin(2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 236 ! a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* & 237 ! (zhin(2)-hginn_d)/2.0) , epsi06) 238 ! ht_i(ji,jj,jl) = zidto(ji,jj) * zhin(2) 239 ! v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 240 ! END DO 241 ! hi_max(ice_cat_bounds(2,1)-1) = zdummy 242 243 ! !rafted ice 244 ! jl = 6 245 ! a_i(ji,jj,jl) = 0.0 246 ! ht_i(ji,jj,jl) = 0.0 247 ! v_i(ji,jj,jl) = 0.0 248 ! END retour a LIMA_MEC 249 !!! 250 251 DO jl = 1, jpl 252 253 !------------- 254 ! Snow depth 255 !------------- 256 ht_s(ji,jj,jl) = zidto(ji,jj) * hninn 257 v_s(ji,jj,jl) = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 258 259 !--------------- 260 ! Ice salinity 261 !--------------- 262 sm_i(ji,jj,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1 263 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 264 265 !---------- 266 ! Ice age 267 !---------- 268 o_i(ji,jj,jl) = zidto(ji,jj) * 1.0 + ( 1.0 - zidto(ji,jj) ) 269 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) 270 271 !------------------------------ 272 ! Sea ice surface temperature 273 !------------------------------ 274 275 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 276 277 !------------------------------------ 278 ! Snow temperature and heat content 279 !------------------------------------ 280 281 DO jk = 1, nlay_s 282 t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 283 ! Snow energy of melting 284 e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 285 ! Change dimensions 286 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 287 ! Multiply by volume, so that heat content in 10^9 Joules 288 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 289 v_s(ji,jj,jl) / nlay_s 290 END DO !jk 291 292 !----------------------------------------------- 293 ! Ice salinities, temperature and heat content 294 !----------------------------------------------- 295 296 DO jk = 1, nlay_i 297 t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 298 s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1 299 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 300 301 ! heat content per unit volume 302 e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 303 ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 304 + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 305 - rcp * ( ztmelts - rtt ) & 306 ) 307 308 ! Correct dimensions to avoid big values 309 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 310 311 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 312 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 313 area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 314 nlay_i 315 END DO ! jk 316 317 END DO ! jl 318 319 ELSE ! on fcor 320 321 !--- Southern hemisphere 322 !---------------------------------------------------------------- 323 324 !----------------------- 325 ! Ice area / thickness 326 !----------------------- 327 328 IF ( jpl .EQ. 1) THEN ! one category 329 330 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 331 a_i(ji,jj,jl) = zidto(ji,jj) * agins_u 332 ht_i(ji,jj,jl) = zidto(ji,jj) * hgins_u 333 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 334 END DO 335 336 ELSE ! several categories 337 338 !level ice 339 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories 340 341 zhis(1) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 342 a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * & 343 (zhis(1)-hgins_u)/2.0) , epsi06 ) 344 ! new line square distribution volume conserving 345 a_i(ji,jj,jl) = zidto(ji,jj) * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) ) 346 ht_i(ji,jj,jl) = zidto(ji,jj) * zhis(1) 347 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 348 349 END DO ! jl 350 351 ENDIF 352 353 !!! 354 ! retour a LIMA_MEC 355 ! !ridged ice 356 ! zdummy = hi_max(ice_cat_bounds(2,1)-1) 357 ! hi_max(ice_cat_bounds(2,1)-1) = 0.0 358 ! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories 359 ! zhis(2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 360 ! a_i(ji,jj,jl) = zidto(ji,jj)*MAX( zgfactors(2) & 361 ! & * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 ) 362 ! ht_i(ji,jj,jl) = zidto(ji,jj) * zhis(2) 363 ! v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 364 ! END DO 365 ! hi_max(ice_cat_bounds(2,1)-1) = zdummy 366 367 ! !rafted ice 368 ! jl = 6 369 ! a_i(ji,jj,jl) = 0.0 370 ! ht_i(ji,jj,jl) = 0.0 371 ! v_i(ji,jj,jl) = 0.0 372 ! END retour a LIMA_MEC 373 !!! 374 375 DO jl = 1, jpl !over thickness categories 376 377 !--------------- 378 ! Snow depth 379 !--------------- 380 381 ht_s(ji,jj,jl) = zidto(ji,jj) * hnins 382 v_s(ji,jj,jl) = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 383 384 !--------------- 385 ! Ice salinity 386 !--------------- 387 388 sm_i(ji,jj,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1 389 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 390 391 !---------- 392 ! Ice age 393 !---------- 394 395 o_i(ji,jj,jl) = zidto(ji,jj) * 1.0 + ( 1.0 - zidto(ji,jj) ) 396 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) 397 398 !------------------------------ 399 ! Sea ice surface temperature 400 !------------------------------ 401 402 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 403 404 !---------------------------------- 405 ! Snow temperature / heat content 406 !---------------------------------- 407 408 DO jk = 1, nlay_s 409 t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 410 ! Snow energy of melting 411 e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 412 ! Change dimensions 413 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 414 ! Multiply by volume, so that heat content in 10^9 Joules 415 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 416 v_s(ji,jj,jl) / nlay_s 417 END DO 418 419 !--------------------------------------------- 420 ! Ice temperature, salinity and heat content 421 !--------------------------------------------- 422 423 DO jk = 1, nlay_i 424 t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 425 s_i(ji,jj,jk,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1 426 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 427 428 ! heat content per unit volume 429 e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 430 ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 431 + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 432 - rcp * ( ztmelts - rtt ) & 433 ) 434 435 ! Correct dimensions to avoid big values 436 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 437 438 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 439 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 440 area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 441 nlay_i 442 END DO !jk 443 444 END DO ! jl 445 446 ENDIF ! on fcor 447 138 zhemis(ji,jj) = 1 ! Northern hemisphere 139 ELSE 140 zhemis(ji,jj) = 2 ! Southern hemisphere 141 ENDIF 448 142 END DO 449 143 END DO 450 451 !-------------------------------------------------------------------- 452 ! 3) Global ice variables for output diagnostics | 453 !-------------------------------------------------------------------- 454 455 fsbbq (:,:) = 0.e0 456 u_ice (:,:) = 0.e0 457 v_ice (:,:) = 0.e0 458 stress1_i(:,:) = 0.0 459 stress2_i(:,:) = 0.0 460 stress12_i(:,:) = 0.0 461 462 !-------------------------------------------------------------------- 463 ! 4) Moments for advection 464 !-------------------------------------------------------------------- 465 466 sxopw (:,:) = 0.e0 467 syopw (:,:) = 0.e0 468 sxxopw(:,:) = 0.e0 469 syyopw(:,:) = 0.e0 470 sxyopw(:,:) = 0.e0 471 472 sxice (:,:,:) = 0.e0 ; sxsn (:,:,:) = 0.e0 ; sxa (:,:,:) = 0.e0 473 syice (:,:,:) = 0.e0 ; sysn (:,:,:) = 0.e0 ; sya (:,:,:) = 0.e0 474 sxxice(:,:,:) = 0.e0 ; sxxsn(:,:,:) = 0.e0 ; sxxa (:,:,:) = 0.e0 475 syyice(:,:,:) = 0.e0 ; syysn(:,:,:) = 0.e0 ; syya (:,:,:) = 0.e0 476 sxyice(:,:,:) = 0.e0 ; sxysn(:,:,:) = 0.e0 ; sxya (:,:,:) = 0.e0 477 478 sxc0 (:,:,:) = 0.e0 ; sxe (:,:,:,:)= 0.e0 479 syc0 (:,:,:) = 0.e0 ; sye (:,:,:,:)= 0.e0 480 sxxc0 (:,:,:) = 0.e0 ; sxxe (:,:,:,:)= 0.e0 481 syyc0 (:,:,:) = 0.e0 ; syye (:,:,:,:)= 0.e0 482 sxyc0 (:,:,:) = 0.e0 ; sxye (:,:,:,:)= 0.e0 483 484 sxsal (:,:,:) = 0.e0 485 sysal (:,:,:) = 0.e0 486 sxxsal (:,:,:) = 0.e0 487 syysal (:,:,:) = 0.e0 488 sxysal (:,:,:) = 0.e0 489 490 !-------------------------------------------------------------------- 491 ! 5) Lateral boundary conditions | 144 ! END MV 2011 new initialization 145 146 !-------------------------------------------------------------------- 147 ! 3) Initialization of sea ice state variables 148 !-------------------------------------------------------------------- 149 150 !----------------------------- 151 ! 3.1) Hemisphere-dependent arrays 152 !----------------------------- 153 ! assign initial thickness, concentration, snow depth and salinity to 154 ! an hemisphere-dependent array 155 zhm_i_ini(1) = hginn ; zhm_i_ini(2) = hgins ! ice thickness 156 zat_i_ini(1) = aginn ; zat_i_ini(2) = agins ! ice concentration 157 zvt_i_ini(:) = zhm_i_ini(:) * zat_i_ini(:) ! ice volume 158 zhm_s_ini(1) = hninn ; zhm_s_ini(2) = hnins ! snow depth 159 zsm_i_ini(1) = sinn ; zsm_i_ini(2) = sins ! bulk ice salinity 160 161 !--------------------------------------------------------------------- 162 ! 3.2) Distribute ice concentration and thickness into the categories 163 !--------------------------------------------------------------------- 164 ! a gaussian distribution for ice concentration is used 165 ! then we check whether the distribution fullfills 166 ! volume and area conservation, positivity and ice categories bounds 167 DO i_hemis = 1, 2 168 169 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 170 171 ! note for the great nemo engineers: 172 ! only very few of the WRITE statements are necessary for the reference version 173 ! they were one day useful, but now i personally doubt of their 174 ! potential for bringing anything useful 175 176 DO i_fill = jpl, 1, -1 177 IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 178 !---------------------------- 179 ! fill the i_fill categories 180 !---------------------------- 181 ! *** 1 category to fill 182 IF ( i_fill .EQ. 1 ) THEN 183 zht_i_ini(1,i_hemis) = zhm_i_ini(i_hemis) 184 za_i_ini(1,i_hemis) = zat_i_ini(i_hemis) 185 zht_i_ini(2:jpl,i_hemis) = 0._wp 186 za_i_ini(2:jpl,i_hemis) = 0._wp 187 ELSE 188 189 ! *** >1 categores to fill 190 !--- Ice thicknesses in the i_fill - 1 first categories 191 DO jl = 1, i_fill - 1 192 zht_i_ini(jl,i_hemis) = 0.5 * ( hi_max(jl) + hi_max(jl-1) ) 193 END DO 194 195 !--- jl0: most likely index where cc will be maximum 196 DO jl = 1, jpl 197 IF ( ( zhm_i_ini(i_hemis) .GT. hi_max(jl-1) ) .AND. & 198 ( zhm_i_ini(i_hemis) .LE. hi_max(jl) ) ) THEN 199 jl0 = jl 200 ENDIF 201 END DO 202 jl0 = MIN(jl0, i_fill) 203 204 !--- Concentrations 205 za_i_ini(jl0,i_hemis) = zat_i_ini(i_hemis) / SQRT(REAL(jpl)) 206 DO jl = 1, i_fill - 1 207 IF ( jl .NE. jl0 ) THEN 208 zsigma = 0.5 * zhm_i_ini(i_hemis) 209 zarg = ( zht_i_ini(jl,i_hemis) - zhm_i_ini(i_hemis) ) / zsigma 210 za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2) 211 ENDIF 212 END DO 213 214 zA = 0. ! sum of the areas in the jpl categories 215 DO jl = 1, i_fill - 1 216 zA = zA + za_i_ini(jl,i_hemis) 217 END DO 218 za_i_ini(i_fill,i_hemis) = zat_i_ini(i_hemis) - zA ! ice conc in the last category 219 IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 220 221 !--- Ice thickness in the last category 222 zV = 0. ! sum of the volumes of the N-1 categories 223 DO jl = 1, i_fill - 1 224 zV = zV + za_i_ini(jl,i_hemis)*zht_i_ini(jl,i_hemis) 225 END DO 226 zht_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 227 IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 228 229 !--- volumes 230 zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zht_i_ini(:,i_hemis) 231 IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 232 233 ENDIF ! i_fill 234 235 !--------------------- 236 ! Compatibility tests 237 !--------------------- 238 ! Test 1: area conservation 239 zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons ) 240 IF ( zconv .LT. 1.0e-6 ) THEN 241 ztest_1 = 1 242 ELSE 243 ! this write is useful 244 IF(lwp) WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 245 ztest_1 = 0 246 ENDIF 247 248 ! Test 2: volume conservation 249 zV_cons = SUM(zv_i_ini(:,i_hemis)) 250 zconv = ABS(zvt_i_ini(i_hemis) - zV_cons) 251 252 IF ( zconv .LT. 1.0e-6 ) THEN 253 ztest_2 = 1 254 ELSE 255 ! this write is useful 256 IF(lwp) WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 257 ' zvt_i_ini = ', zvt_i_ini(i_hemis) 258 ztest_2 = 0 259 ENDIF 260 261 ! Test 3: thickness of the last category is in-bounds ? 262 IF ( zht_i_ini(i_fill, i_hemis) .GT. hi_max(i_fill-1) ) THEN 263 ztest_3 = 1 264 ELSE 265 ! this write is useful 266 IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zht_i_ini(i_fill,i_hemis) = ', & 267 zht_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 268 ztest_3 = 0 269 ENDIF 270 271 ! Test 4: positivity of ice concentrations 272 ztest_4 = 1 273 DO jl = 1, jpl 274 IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 275 ! this write is useful 276 IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis) 277 ztest_4 = 0 278 ENDIF 279 END DO 280 281 ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 282 283 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 284 285 END DO ! i_fill 286 287 IF(lwp) THEN 288 WRITE(numout,*), ' ztests : ', ztests 289 IF ( ztests .NE. 4 ) THEN 290 WRITE(numout,*) 291 WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 292 WRITE(numout,*), ' !!!! RED ALERT !!! ' 293 WRITE(numout,*), ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!' 294 WRITE(numout,*), ' !!!! Something is wrong in the LIM3 initialization procedure ' 295 WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 296 WRITE(numout,*) 297 WRITE(numout,*), ' *** ztests is not equal to 4 ' 298 WRITE(numout,*), ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 299 WRITE(numout,*), ' zat_i_ini : ', zat_i_ini(i_hemis) 300 WRITE(numout,*), ' zhm_i_ini : ', zhm_i_ini(i_hemis) 301 ENDIF ! ztests .NE. 4 302 ENDIF 303 304 END DO ! i_hemis 305 306 !--------------------------------------------------------------------- 307 ! 3.3) Space-dependent arrays for ice state variables 308 !--------------------------------------------------------------------- 309 310 ! Ice concentration, thickness and volume, snow depth, ice 311 ! salinity, ice age, surface temperature 312 DO jl = 1, jpl ! loop over categories 313 DO jj = 1, jpj 314 DO ji = 1, jpi 315 a_i(ji,jj,jl) = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj)) ! concentration 316 ht_i(ji,jj,jl) = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj)) ! ice thickness 317 ht_s(ji,jj,jl) = zidto(ji,jj) * zhm_s_ini(zhemis(ji,jj)) ! snow depth 318 sm_i(ji,jj,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) ! salinity 319 o_i(ji,jj,jl) = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age 320 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * t_bo(ji,jj) ! surf temp 321 322 ! ice volume, snow volume, salt content, age content 323 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume 324 v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume 325 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 326 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content 327 END DO ! ji 328 END DO ! jj 329 END DO ! jl 330 331 ! Snow temperature and heat content 332 DO jk = 1, nlay_s 333 DO jl = 1, jpl ! loop over categories 334 DO jj = 1, jpj 335 DO ji = 1, jpi 336 t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * rtt 337 ! Snow energy of melting 338 e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 339 ! Change dimensions 340 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 341 ! Multiply by volume, so that heat content in 10^9 Joules 342 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 343 END DO ! ji 344 END DO ! jj 345 END DO ! jl 346 END DO ! jk 347 348 ! Ice salinity, temperature and heat content 349 DO jk = 1, nlay_i 350 DO jl = 1, jpl ! loop over categories 351 DO jj = 1, jpj 352 DO ji = 1, jpi 353 t_i(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1._wp - zidto(ji,jj) ) * rtt 354 s_i(ji,jj,jk,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min 355 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 356 357 ! heat content per unit volume 358 e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 359 + lfus * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 360 - rcp * ( ztmelts - rtt ) ) 361 362 ! Correct dimensions to avoid big values 363 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 364 365 ! Mutliply by ice volume, and divide by number of layers 366 ! to get heat content in 10^9 J 367 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i 368 END DO ! ji 369 END DO ! jj 370 END DO ! jl 371 END DO ! jk 372 373 !-------------------------------------------------------------------- 374 ! 4) Global ice variables for output diagnostics | 375 !-------------------------------------------------------------------- 376 fsbbq (:,:) = 0._wp 377 u_ice (:,:) = 0._wp 378 v_ice (:,:) = 0._wp 379 stress1_i(:,:) = 0._wp 380 stress2_i(:,:) = 0._wp 381 stress12_i(:,:) = 0._wp 382 383 # if defined key_coupled 384 albege(:,:) = 0.8 * tms(:,:) 385 # endif 386 387 !-------------------------------------------------------------------- 388 ! 5) Moments for advection 389 !-------------------------------------------------------------------- 390 391 sxopw (:,:) = 0._wp 392 syopw (:,:) = 0._wp 393 sxxopw(:,:) = 0._wp 394 syyopw(:,:) = 0._wp 395 sxyopw(:,:) = 0._wp 396 397 sxice (:,:,:) = 0._wp ; sxsn (:,:,:) = 0._wp ; sxa (:,:,:) = 0._wp 398 syice (:,:,:) = 0._wp ; sysn (:,:,:) = 0._wp ; sya (:,:,:) = 0._wp 399 sxxice(:,:,:) = 0._wp ; sxxsn(:,:,:) = 0._wp ; sxxa (:,:,:) = 0._wp 400 syyice(:,:,:) = 0._wp ; syysn(:,:,:) = 0._wp ; syya (:,:,:) = 0._wp 401 sxyice(:,:,:) = 0._wp ; sxysn(:,:,:) = 0._wp ; sxya (:,:,:) = 0._wp 402 403 sxc0 (:,:,:) = 0._wp ; sxe (:,:,:,:)= 0._wp 404 syc0 (:,:,:) = 0._wp ; sye (:,:,:,:)= 0._wp 405 sxxc0 (:,:,:) = 0._wp ; sxxe (:,:,:,:)= 0._wp 406 syyc0 (:,:,:) = 0._wp ; syye (:,:,:,:)= 0._wp 407 sxyc0 (:,:,:) = 0._wp ; sxye (:,:,:,:)= 0._wp 408 409 sxsal (:,:,:) = 0._wp 410 sysal (:,:,:) = 0._wp 411 sxxsal (:,:,:) = 0._wp 412 syysal (:,:,:) = 0._wp 413 sxysal (:,:,:) = 0._wp 414 415 !-------------------------------------------------------------------- 416 ! 6) Lateral boundary conditions | 492 417 !-------------------------------------------------------------------- 493 418 494 419 DO jl = 1, jpl 420 495 421 CALL lbc_lnk( a_i(:,:,jl) , 'T', 1. ) 496 422 CALL lbc_lnk( v_i(:,:,jl) , 'T', 1. ) … … 498 424 CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. ) 499 425 CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. ) 500 ! 426 501 427 CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. ) 502 428 CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. ) … … 515 441 a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl) 516 442 END DO 443 444 at_i (:,:) = 0.0_wp 445 DO jl = 1, jpl 446 at_i (:,:) = at_i (:,:) + a_i (:,:,jl) 447 END DO 517 448 518 449 CALL lbc_lnk( at_i , 'T', 1. ) … … 521 452 CALL lbc_lnk( fsbbq , 'T', 1. ) 522 453 ! 523 CALL wrk_dealloc( jpm, zgfactorn, zgfactors, zhin, zhis ) 454 !-------------------------------------------------------------------- 455 ! 6) ???? | 456 !-------------------------------------------------------------------- 457 tn_ice (:,:,:) = t_su (:,:,:) 458 524 459 CALL wrk_dealloc( jpi, jpj, zidto ) 525 ! 460 CALL wrk_dealloc( jpi, jpj, zhemis ) 461 CALL wrk_dealloc( jpl, 2, zht_i_ini, za_i_ini, zv_i_ini ) 462 CALL wrk_dealloc( 2, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 463 526 464 END SUBROUTINE lim_istate 527 528 465 529 466 SUBROUTINE lim_istate_init … … 533 470 !! ** Purpose : Definition of initial state of the ice 534 471 !! 535 !! ** Method : Read the namiceini namelist and check the parameter 536 !! values called at the first timestep (nit000) 537 !! 538 !! ** input : namelist namiceini 472 !! ** Method : Read the namiceini namelist and check the parameter 473 !! values called at the first timestep (nit000) 474 !! 475 !! ** input : 476 !! Namelist namiceini 477 !! 478 !! history : 479 !! 8.5 ! 03-08 (C. Ethe) original code 480 !! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization 539 481 !!----------------------------------------------------------------------------- 540 NAMELIST/namiceini/ ttest, hninn, hginn_u, aginn_u, hginn_d, aginn_d, hnins, & 541 & hgins_u, agins_u, hgins_d, agins_d, sinn, sins 482 NAMELIST/namiceini/ ttest, hninn, hnins, hginn, hgins, aginn, agins, sinn, sins 542 483 !!----------------------------------------------------------------------------- 543 ! 544 REWIND ( numnam_ice ) ! Read Namelist namiceini 484 485 ! Define the initial parameters 486 ! ------------------------- 487 488 ! Read Namelist namiceini 489 REWIND ( numnam_ice ) 545 490 READ ( numnam_ice , namiceini ) 546 ! 547 IF(lwp) THEN ! control print 491 IF(lwp) THEN 548 492 WRITE(numout,*) 549 493 WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' … … 551 495 WRITE(numout,*) ' threshold water temp. for initial sea-ice ttest = ', ttest 552 496 WRITE(numout,*) ' initial snow thickness in the north hninn = ', hninn 553 WRITE(numout,*) ' initial undef ice thickness in the north hginn_u = ', hginn_u554 WRITE(numout,*) ' initial undef ice concentr. in the north aginn_u = ', aginn_u555 WRITE(numout,*) ' initial def ice thickness in the north hginn_d = ', hginn_d556 WRITE(numout,*) ' initial def ice concentr. in the north aginn_d = ', aginn_d557 497 WRITE(numout,*) ' initial snow thickness in the south hnins = ', hnins 558 WRITE(numout,*) ' initial undef ice thickness in the north hgins_u = ', hgins_u559 WRITE(numout,*) ' initial undef ice concentr. in the north agins_u = ', agins_u560 WRITE(numout,*) ' initial def ice thickness in the north hgins_d = ', hgins_d561 WRITE(numout,*) ' initial def ice concentr. in the north agins_d = ', agins_d562 WRITE(numout,*) ' initial ice salinity in the northsinn = ', sinn563 WRITE(numout,*) ' initial ice salinity in the southsins = ', sins498 WRITE(numout,*) ' initial ice thickness in the north hginn = ', hginn 499 WRITE(numout,*) ' initial ice thickness in the south hgins = ', hgins 500 WRITE(numout,*) ' initial ice concentr. in the north aginn = ', aginn 501 WRITE(numout,*) ' initial ice concentr. in the north agins = ', agins 502 WRITE(numout,*) ' initial ice salinity in the north sinn = ', sinn 503 WRITE(numout,*) ' initial ice salinity in the south sins = ', sins 564 504 ENDIF 565 ! 505 566 506 END SUBROUTINE lim_istate_init 567 507 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r4045 r4099 416 416 417 417 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 ) 418 deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 419 !!gm faster to replace the line above with simply: 420 !! deltat(ji,jj) = MAX( delta, creepl ) 421 !!gm end 422 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 423 425 !-Calculate stress tensor components zs1 and zs2 424 426 !-at centre of grid cells (see section 3.5 of CICE user's guide). … … 741 743 & - e1v( ji , jj-1 ) * u_ice2(ji ,jj-1) ) / area(ji,jj) 742 744 743 deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) & 744 & + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 745 & ) + creepl 746 745 ! deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) & 746 ! & + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 747 ! & ) + creepl 748 ! MV rewriting 749 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 ) 750 deltat(ji,jj) = delta + creepl 751 ! END MV 752 747 753 ENDIF ! zdummy 748 754 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
r4045 r4099 371 371 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 372 372 t_su(:,:,jl) = z2d(:,:) 373 tn_ice (:,:,:) = t_su (:,:,:) 373 374 END DO 374 375 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r4045 r4099 102 102 INTEGER, INTENT(in) :: kt ! number of iteration 103 103 ! 104 INTEGER :: ji, jj ! dummy loop indices104 INTEGER :: ji, jj, jl ! dummy loop indices 105 105 INTEGER :: ierr, ifvt, i1mfr, idfr ! local integer 106 106 INTEGER :: iflt, ial , iadv , ifral, ifrdv ! - - … … 145 145 146 146 ! computation the solar flux at ocean surface 147 zfcm1 = pfrld(ji,jj) * qsr(ji,jj) + & 148 & ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 147 IF (lk_cpl) THEN ! be carfeful: not being tested yet 148 ! original line 149 !zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) 150 ! new line to include solar penetration (not tested) 151 zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 152 DO jl = 1, jpl 153 zfcm1 = zfcm1 - qsr_ice(ji,jj,jl) * a_i(ji,jj,jl) 154 END DO 155 ELSE 156 zfcm1 = pfrld(ji,jj) * qsr(ji,jj) + & 157 & ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 158 ENDIF 149 159 ! fstric Solar flux transmitted trough the ice 150 160 ! qsr Net short wave heat flux on free ocean … … 227 237 228 238 ! computing freshwater exchanges at the ice/ocean interface 229 zemp = emp(ji,jj) * ( 1.0 - at_i(ji,jj) ) & ! evaporation over oceanic fraction 230 & - tprecip(ji,jj) * at_i(ji,jj) & ! all precipitation reach the ocean 231 & + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) ) & ! except solid precip intercepted by sea-ice 232 & - fmmec(ji,jj) ! snow falling when ridging 239 IF (lk_cpl) THEN 240 zemp = - emp_tot(ji,jj) + emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) ) & ! 241 & - rdm_snw(ji,jj) / rdt_ice 242 ELSE 243 zemp = emp(ji,jj) * ( 1.0 - at_i(ji,jj) ) & ! evaporation over oceanic fraction 244 & - tprecip(ji,jj) * at_i(ji,jj) & ! all precipitation reach the ocean 245 & + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) ) & ! except solid precip intercepted by sea-ice 246 & - fmmec(ji,jj) ! snow falling when ridging 247 ENDIF 233 248 234 249 ! mass flux at the ocean/ice interface (sea ice fraction) -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4045 r4099 84 84 REAL(wp) :: zdhnm, zhnnew, zhisn, zihic, zzc ! 85 85 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 86 REAL(wp) :: zds ! increment of bottom ice salinity87 86 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads 88 87 REAL(wp) :: zsm_snowice ! snow-ice salinity … … 429 428 zfracs = zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & 430 429 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 431 z ds = zfracs * sss_m(ii,ij) - s_i_new(ji)430 zfracs = MIN( 0.5 , zfracs ) 432 431 s_i_new(ji) = zfracs * sss_m(ii,ij) 433 432 ENDIF ! fc_bo_i … … 438 437 DO ji = kideb, kiut 439 438 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN 440 ! New ice salinity must not exceed 15psu439 ! New ice salinity must not exceed 20 psu 441 440 s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 442 441 ! Metling point in K -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r4072 r4099 26 26 USE limthd 27 27 USE limsbc 28 USE limdia29 28 USE limdiahsb 30 29 USE limwri … … 402 401 ! 2.11) Ice salinity 403 402 !--------------------- 404 !clem@bug: smv_i should be updated too: smv_i(:,:,:) = smv_i(:,:,:) + sm_i(:,:,:) * ( v_i(:,:,:) - zviold(:,:,:) ) 403 ! clem correct bug on smv_i 404 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 405 405 406 IF ( num_sal == 2 ) THEN ! general case 406 407 407 DO jl = 1, jpl 408 408 !DO jk = 1, nlay_i … … 410 410 DO ji = 1, jpi 411 411 ! salinity stays in bounds 412 smv_i(ji,jj,jl) = MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) ) 413 i_ice_switch = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl) + epsi20)) 414 smv_i(ji,jj,jl) = i_ice_switch*smv_i(ji,jj,jl) + 0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl) 412 !clem smv_i(ji,jj,jl) = MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) ) 413 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 414 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20 ) ) 415 smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 415 416 END DO ! ji 416 417 END DO ! jj 417 418 !END DO !jk 418 419 END DO !jl 419 420 420 ENDIF 421 421 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r4072 r4099 26 26 USE limthd 27 27 USE limsbc 28 USE limdia29 28 USE limdiahsb 30 29 USE limwri … … 305 304 zesum = zesum + e_i(ji,jj,jk,jl) 306 305 END DO 307 IF (ind_im .LT. nlay_i ) THEN308 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) )309 ENDIF306 !clem IF (ind_im .LT. nlay_i ) THEN 307 !clem smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) ) 308 !clem ENDIF 310 309 ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) 311 310 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) … … 532 531 ! 2.11) Ice salinity 533 532 !--------------------- 534 !clem@bug: smv_i should be updated too: smv_i(:,:,:) = smv_i(:,:,:) + sm_i(:,:,:) * ( v_i(:,:,:) - zviold(:,:,:) ) 533 ! clem correct bug on smv_i 534 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 535 535 536 IF ( num_sal == 2 ) THEN ! general case 536 537 DO jl = 1, jpl … … 539 540 DO ji = 1, jpi 540 541 ! salinity stays in bounds 541 smv_i(ji,jj,jl) = MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) ) 542 i_ice_switch = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl) + epsi20)) 543 smv_i(ji,jj,jl) = i_ice_switch*smv_i(ji,jj,jl) + 0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl) 542 !clem smv_i(ji,jj,jl) = MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) ) 543 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 544 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20 ) ) 545 smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 544 546 END DO ! ji 545 547 END DO ! jj -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
r3625 r4099 30 30 31 31 PUBLIC sbc_ice_alloc ! called in iceini(_2).F90 32 33 CHARACTER (len=8), PUBLIC :: cn_iceflx = 'none' !: Flux handling over ice categories 34 LOGICAL, PUBLIC :: ln_iceflx_ave = .FALSE. ! Average heat fluxes over all ice categories 35 LOGICAL, PUBLIC :: ln_iceflx_linear = .FALSE. ! Redistribute mean heat fluxes over all ice categories, using ice temperature and albedo 32 36 33 37 # if defined key_lim2 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
r3914 r4099 471 471 ! Coupled case: since cloud cover is not received from atmosphere 472 472 ! ===> defined as constant value -> definition done in sbc_cpl_init 473 fr1_i0(:,:) = 0.18474 fr2_i0(:,:) = 0.82473 IF ( ALLOCATED (fr1_i0)) fr1_i0 (:,:) = 0.18 474 IF ( ALLOCATED (fr2_i0)) fr2_i0 (:,:) = 0.82 475 475 ! ! ------------------------- ! 476 476 ! ! 10m wind module ! … … 931 931 CALL wrk_alloc( jpi,jpj, ztx, zty ) 932 932 933 IF( srcv(jpr_itx1)%laction ) THEN ; itx = jpr_itx1 933 !AC Pour eviter un stress nul sur la glace dans le cas mixed oce-ice 934 IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN ; itx = jpr_itx1 934 935 ELSE ; itx = jpr_otx1 935 936 ENDIF … … 938 939 IF( nrcvinfo(itx) == OASIS_Rcv ) THEN 939 940 940 ! ! ======================= ! 941 IF( srcv(jpr_itx1)%laction ) THEN ! ice stress received ! 942 ! ! ======================= ! 941 ! ! ======================= ! 942 !AC Pour eviter un stress nul sur la glace dans le cas mixes oce-ice 943 IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN ! ice stress received ! 944 ! ! ======================= ! 943 945 ! 944 946 IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN ! 2 components on the sphere -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90
r3625 r4099 99 99 100 100 fr_i(:,:) = tfreez( sss_m ) * tmask(:,:,1) ! sea surface freezing temperature [Celcius] 101 #if defined key_coupled 101 102 ! OM : probleme. a_i pas defini dans les cas lim3 et cice 103 #if defined key_coupled && defined key_lim2 102 104 a_i(:,:,1) = fr_i(:,:) 103 105 #endif -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r4042 r4099 32 32 USE sbcblk_core ! Surface boundary condition: CORE bulk 33 33 USE sbcblk_clio ! Surface boundary condition: CLIO bulk 34 USE sbccpl ! Surface boundary condition: coupled interface 34 35 USE albedo ! ocean & ice albedo 35 36 … … 42 43 USE limitd_me ! Mechanics on ice thickness distribution 43 44 USE limsbc ! sea surface boundary condition 44 USE limdia ! Ice diagnostics45 45 USE limdiahsb ! Ice budget diagnostics 46 46 USE limwri ! Ice outputs … … 77 77 !!---------------------------------------------------------------------- 78 78 CONTAINS 79 80 FUNCTION fice_cell_ave ( ptab) 81 !!-------------------------------------------------------------------------- 82 !! * Compute average over categories, for grid cell (ice covered and free ocean) 83 !!-------------------------------------------------------------------------- 84 REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave 85 REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab 86 INTEGER :: jl ! Dummy loop index 87 88 fice_cell_ave (:,:) = 0.0_wp 89 90 DO jl = 1, jpl 91 fice_cell_ave (:,:) = fice_cell_ave (:,:) & 92 & + a_i (:,:,jl) * ptab (:,:,jl) 93 END DO 94 95 END FUNCTION fice_cell_ave 96 97 FUNCTION fice_ice_ave ( ptab) 98 !!-------------------------------------------------------------------------- 99 !! * Compute average over categories, for ice covered part of grid cell 100 !!-------------------------------------------------------------------------- 101 REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave 102 REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab 103 104 fice_ice_ave (:,:) = 0.0_wp 105 WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:) 106 107 END FUNCTION fice_ice_ave 108 109 !!====================================================================== 79 110 80 111 SUBROUTINE sbc_ice_lim( kt, kblk ) … … 104 135 REAL(wp) :: zcoef ! local scalar 105 136 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_ice_os, zalb_ice_cs ! albedo of the ice under overcast/clear sky 137 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_ice ! mean albedo of ice (for coupled) 138 139 REAL(wp), POINTER, DIMENSION(:,:) :: zalb_ice_all ! Mean albedo over all categories 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztem_ice_all ! Mean temperature over all categories 141 142 REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_ice_all ! Mean solar heat flux over all categories 143 REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_ice_all ! Mean non solar heat flux over all categories 144 REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_ice_all ! Mean latent heat flux over all categories 145 REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all ! Mean d(qns)/dT over all categories 146 REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all ! Mean d(qla)/dT over all categories 106 147 !!---------------------------------------------------------------------- 107 148 149 !- O.M. : why do we allocate all these arrays even when MOD( kt-1, nn_fsbc ) /= 0 ????? 150 108 151 IF( nn_timing == 1 ) CALL timing_start('sbc_ice_lim') 109 152 110 153 CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 154 155 IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 156 CALL wrk_alloc( jpi,jpj,jpl, zalb_ice) 157 END IF 158 IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 159 CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 160 ENDIF 161 111 162 112 163 IF( kt == nit000 ) THEN … … 139 190 t_su(:,:,jl) = t_su(:,:,jl) + rt0 * ( 1. - tmask(:,:,1) ) 140 191 END DO 192 193 IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) + zalb_ice_os (:,:,:) ) 194 195 IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 196 ! 197 ! Compute mean albedo and temperature 198 zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) ) 199 ztem_ice_all (:,:) = fice_ice_ave ( tn_ice (:,:,:) ) 200 ! 201 ENDIF 141 202 ! Bulk formulea - provides the following fields: 142 203 ! utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2] … … 161 222 & tprecip , sprecip , & 162 223 & fr1_i0 , fr2_i0 , cp_ice_msh, jpl ) 224 ! 225 CASE ( 5 ) 226 zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) + zalb_ice_os (:,:,:) ) 227 228 CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) 229 230 CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=tn_ice ) 231 232 ! Latent heat flux is forced to 0 in coupled : 233 ! it is included in qns (non-solar heat flux) 234 qla_ice (:,:,:) = 0.0e0_wp 235 dqla_ice (:,:,:) = 0.0e0_wp 236 ! 163 237 END SELECT 238 239 ! Average over all categories 240 IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 241 242 z_qns_ice_all (:,:) = fice_ice_ave ( qns_ice (:,:,:) ) 243 z_qsr_ice_all (:,:) = fice_ice_ave ( qsr_ice (:,:,:) ) 244 z_dqns_ice_all (:,:) = fice_ice_ave ( dqns_ice (:,:,:) ) 245 z_qla_ice_all (:,:) = fice_ice_ave ( qla_ice (:,:,:) ) 246 z_dqla_ice_all (:,:) = fice_ice_ave ( dqla_ice (:,:,:) ) 247 248 DO jl = 1, jpl 249 dqns_ice (:,:,jl) = z_dqns_ice_all (:,:) 250 dqla_ice (:,:,jl) = z_dqla_ice_all (:,:) 251 END DO 252 ! 253 IF ( ln_iceflx_ave ) THEN 254 DO jl = 1, jpl 255 qns_ice (:,:,jl) = z_qns_ice_all (:,:) 256 qsr_ice (:,:,jl) = z_qsr_ice_all (:,:) 257 qla_ice (:,:,jl) = z_qla_ice_all (:,:) 258 END DO 259 END IF 260 ! 261 IF ( ln_iceflx_linear ) THEN 262 DO jl = 1, jpl 263 qns_ice (:,:,jl) = z_qns_ice_all(:,:) + z_dqns_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:)) 264 qla_ice (:,:,jl) = z_qla_ice_all(:,:) + z_dqla_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:)) 265 qsr_ice (:,:,jl) = (1.0e0_wp-zalb_ice(:,:,jl)) / (1.0e0_wp-zalb_ice_all(:,:)) * z_qsr_ice_all(:,:) 266 END DO 267 END IF 268 END IF 164 269 165 270 ! !----------------------! … … 264 369 ! 265 370 ! ! Diagnostics and outputs 266 IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp ) & 267 & CALL lim_dia 268 IF (ln_limdiahsb) CALL lim_diahsb 371 IF (ln_limdiaout) CALL lim_diahsb 372 !clem # if ! defined key_iomput 269 373 CALL lim_wri( 1 ) ! Ice outputs 374 !clem # endif 270 375 IF( kt == nit000 ) CALL iom_close( numrir ) ! clem: close input ice restart file 271 376 IF( lrst_ice ) CALL lim_rst_write( kt ) ! Ice restart file … … 287 392 ! 288 393 CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 394 IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 395 CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice) 396 END IF 397 IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 398 CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 399 ENDIF 289 400 ! 290 401 IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_lim') -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r4054 r4099 84 84 NAMELIST/namsbc/ nn_fsbc , ln_ana , ln_flx, ln_blk_clio, ln_blk_core, ln_cpl, & 85 85 & ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc , ln_rnf, & 86 & ln_ssr , nn_fwb , ln_cdgw , ln_wave , ln_sdw 86 & ln_ssr , nn_fwb , ln_cdgw , ln_wave , ln_sdw, cn_iceflx 87 87 !!---------------------------------------------------------------------- 88 88 … … 117 117 WRITE(numout,*) ' MFS bulk formulation ln_blk_mfs = ', ln_blk_mfs 118 118 WRITE(numout,*) ' coupled formulation (T if key_sbc_cpl) ln_cpl = ', ln_cpl 119 WRITE(numout,*) ' Flux handling over ice categories cn_iceflx = ', TRIM (cn_iceflx) 119 120 WRITE(numout,*) ' Misc. options of sbc : ' 120 121 WRITE(numout,*) ' Patm gradient added in ocean & ice Eqs. ln_apr_dyn = ', ln_apr_dyn … … 128 129 ENDIF 129 130 131 ! Flux handling over ice categories 132 SELECT CASE ( TRIM (cn_iceflx)) 133 CASE ('ave') 134 ln_iceflx_ave = .TRUE. 135 ln_iceflx_linear = .FALSE. 136 CASE ('linear') 137 ln_iceflx_ave = .FALSE. 138 ln_iceflx_linear = .TRUE. 139 CASE default 140 ln_iceflx_ave = .FALSE. 141 ln_iceflx_linear = .FALSE. 142 END SELECT 143 IF(lwp) WRITE(numout,*) ' Fluxes averaged over all ice categories ln_iceflx_ave = ', ln_iceflx_ave 144 IF(lwp) WRITE(numout,*) ' Fluxes distributed linearly over ice categories ln_iceflx_linear = ', ln_iceflx_linear 145 ! 130 146 ! ! allocate sbc arrays 131 147 IF( sbc_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' ) … … 166 182 IF( ( nn_ice == 3 .OR. nn_ice == 4 ) .AND. nn_ice_embd == 0 ) & 167 183 & CALL ctl_stop( 'LIM3 and CICE sea-ice models require nn_ice_embd = 1 or 2' ) 184 185 IF( ln_iceflx_ave .AND. ln_iceflx_linear ) & 186 & CALL ctl_stop( ' ln_iceflx_ave and ln_iceflx_linear options are not compatible' ) 187 188 IF( ( nn_ice ==3 .AND. lk_cpl) .AND. .NOT. ( ln_iceflx_ave .OR. ln_iceflx_linear ) ) & 189 & CALL ctl_stop( ' With lim3 coupled, either ln_iceflx_ave or ln_iceflx_linear must be set to .TRUE.' ) 168 190 169 191 IF( ln_dm2dc ) nday_qsr = -1 ! initialisation flag … … 295 317 CASE( 2 ) ; CALL sbc_ice_lim_2( kt, nsbc ) ! LIM-2 ice model 296 318 CASE( 3 ) ; CALL sbc_ice_lim ( kt, nsbc ) ! LIM-3 ice model 319 !is it useful? 320 !CASE( 4 ) ; CALL sbc_ice_cice ( kt, nsbc ) ! CICE ice model 297 321 END SELECT 298 322
Note: See TracChangeset
for help on using the changeset viewer.
