- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/ICE/iceistate.F90
r11229 r13463 22 22 USE eosbn2 ! equation of state 23 23 USE domvvl ! Variable volume 24 USE ice ! sea-ice variables 25 USE icevar ! ice_var_salprof 24 USE ice ! sea-ice: variables 25 USE ice1D ! sea-ice: thermodynamics variables 26 USE icetab ! sea-ice: 1D <==> 2D transformation 27 USE icevar ! sea-ice: operations 26 28 ! 27 29 USE in_out_manager ! I/O manager … … 31 33 USE fldread ! read input fields 32 34 35 # if defined key_agrif 36 USE agrif_oce 37 USE agrif_ice 38 USE agrif_ice_interp 39 # endif 40 33 41 IMPLICIT NONE 34 42 PRIVATE … … 36 44 PUBLIC ice_istate ! called by icestp.F90 37 45 PUBLIC ice_istate_init ! called by icestp.F90 38 39 INTEGER , PARAMETER :: jpfldi = 6 ! maximum number of files to read40 INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m) at T-point41 INTEGER , PARAMETER :: jp_hts = 2 ! index of snow thicknes (m) at T-point42 INTEGER , PARAMETER :: jp_ati = 3 ! index of ice fraction (%) at T-point43 INTEGER , PARAMETER :: jp_tsu = 4 ! index of ice surface temp (K) at T-point44 INTEGER , PARAMETER :: jp_tmi = 5 ! index of ice temp at T-point45 INTEGER , PARAMETER :: jp_smi = 6 ! index of ice sali at T-point46 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read)47 46 ! 48 47 ! !! ** namelist (namini) ** 49 LOGICAL :: ln_iceini ! initialization or not 50 LOGICAL :: ln_iceini_file ! Ice initialization state from 2D netcdf file 51 REAL(wp) :: rn_thres_sst ! threshold water temperature for initial sea ice 52 REAL(wp) :: rn_hts_ini_n ! initial snow thickness in the north 53 REAL(wp) :: rn_hts_ini_s ! initial snow thickness in the south 54 REAL(wp) :: rn_hti_ini_n ! initial ice thickness in the north 55 REAL(wp) :: rn_hti_ini_s ! initial ice thickness in the south 56 REAL(wp) :: rn_ati_ini_n ! initial leads area in the north 57 REAL(wp) :: rn_ati_ini_s ! initial leads area in the south 58 REAL(wp) :: rn_smi_ini_n ! initial salinity 59 REAL(wp) :: rn_smi_ini_s ! initial salinity 60 REAL(wp) :: rn_tmi_ini_n ! initial temperature 61 REAL(wp) :: rn_tmi_ini_s ! initial temperature 62 48 LOGICAL, PUBLIC :: ln_iceini !: Ice initialization or not 49 LOGICAL, PUBLIC :: ln_iceini_file !: Ice initialization from 2D netcdf file 50 REAL(wp) :: rn_thres_sst 51 REAL(wp) :: rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n 52 REAL(wp) :: rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s 53 REAL(wp) :: rn_apd_ini_n, rn_hpd_ini_n 54 REAL(wp) :: rn_apd_ini_s, rn_hpd_ini_s 55 ! 56 ! ! if ln_iceini_file = T 57 INTEGER , PARAMETER :: jpfldi = 9 ! maximum number of files to read 58 INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m) 59 INTEGER , PARAMETER :: jp_hts = 2 ! index of snw thickness (m) 60 INTEGER , PARAMETER :: jp_ati = 3 ! index of ice fraction (-) 61 INTEGER , PARAMETER :: jp_smi = 4 ! index of ice salinity (g/kg) 62 INTEGER , PARAMETER :: jp_tmi = 5 ! index of ice temperature (K) 63 INTEGER , PARAMETER :: jp_tsu = 6 ! index of ice surface temp (K) 64 INTEGER , PARAMETER :: jp_tms = 7 ! index of snw temperature (K) 65 INTEGER , PARAMETER :: jp_apd = 8 ! index of pnd fraction (-) 66 INTEGER , PARAMETER :: jp_hpd = 9 ! index of pnd depth (m) 67 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read) 68 69 !! * Substitutions 70 # include "do_loop_substitute.h90" 63 71 !!---------------------------------------------------------------------- 64 72 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 68 76 CONTAINS 69 77 70 SUBROUTINE ice_istate 78 SUBROUTINE ice_istate( kt, Kbb, Kmm, Kaa ) 71 79 !!------------------------------------------------------------------- 72 80 !! *** ROUTINE ice_istate *** … … 87 95 !! 88 96 !! ** Notes : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even 89 !! where there is no ice (clem: I do not know why, is it mandatory?)97 !! where there is no ice 90 98 !!-------------------------------------------------------------------- 99 INTEGER, INTENT(in) :: kt ! time step 100 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 101 ! 91 102 INTEGER :: ji, jj, jk, jl ! dummy loop indices 92 INTEGER :: i_hemis, i_fill, jl0 ! local integers 93 REAL(wp) :: ztmelts, zdh 94 REAL(wp) :: zarg, zV, zconv, zdv, zfac 103 REAL(wp) :: ztmelts 95 104 INTEGER , DIMENSION(4) :: itest 96 105 REAL(wp), DIMENSION(jpi,jpj) :: z2d 97 106 REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator 98 REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file 99 REAL(wp), DIMENSION(jpi,jpj) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 100 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini , za_i_ini !data by cattegories to fill 107 REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file 108 REAL(wp), DIMENSION(jpi,jpj) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 109 REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini !data from namelist or nc file 110 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !locak arrays 111 !! 112 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d 101 113 !-------------------------------------------------------------------- 102 114 … … 105 117 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 106 118 107 !-------------------------------------------------------------------- 108 ! 1) Set surface and bottom temperatures to initial values 109 !-------------------------------------------------------------------- 110 ! 111 ! init surface temperature 119 !--------------------------- 120 ! 1) 1st init. of the fields 121 !--------------------------- 122 ! 123 ! basal temperature (considered at freezing point) [Kelvin] 124 CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 125 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 126 ! 127 ! surface temperature and conductivity 112 128 DO jl = 1, jpl 113 129 t_su (:,:,jl) = rt0 * tmask(:,:,1) ! temp at the surface … … 115 131 END DO 116 132 ! 117 ! init basal temperature (considered at freezing point) [Kelvin] 118 CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 119 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 120 121 IF( ln_iceini ) THEN 122 !----------------------------------------------------------- 123 ! 2) Compute or read sea ice variables ===> single category 124 !----------------------------------------------------------- 125 ! 126 ! !---------------! 127 IF( ln_iceini_file )THEN ! Read a file ! 128 ! !---------------! 129 ! 130 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) 131 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) 132 zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) 133 zts_u_ini(:,:) = si(jp_tsu)%fnow(:,:,1) 134 ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) 135 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) 136 ! 137 WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 138 ELSEWHERE ; zswitch(:,:) = 0._wp 139 END WHERE 140 zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) 141 ! 142 ! !---------------! 143 ELSE ! Read namelist ! 144 ! !---------------! 145 ! no ice if sst <= t-freez + ttest 146 WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 147 ELSEWHERE ; zswitch(:,:) = tmask(:,:,1) 148 END WHERE 149 ! 150 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 151 WHERE( ff_t(:,:) >= 0._wp ) 152 zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:) 153 zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 154 zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 155 zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 156 zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 157 ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 158 ELSEWHERE 159 zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 160 zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 161 zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 162 zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 163 zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 164 ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 165 END WHERE 166 zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) 167 ! 168 ENDIF 169 170 !------------------------------------------------------------------ 171 ! 3) Distribute ice concentration and thickness into the categories 172 !------------------------------------------------------------------ 173 ! a gaussian distribution for ice concentration is used 174 ! then we check whether the distribution fullfills 175 ! volume and area conservation, positivity and ice categories bounds 176 177 IF( jpl == 1 ) THEN 178 ! 179 zh_i_ini(:,:,1) = zht_i_ini(:,:) 180 za_i_ini(:,:,1) = zat_i_ini(:,:) 181 ! 182 ELSE 183 zh_i_ini(:,:,:) = 0._wp 184 za_i_ini(:,:,:) = 0._wp 185 ! 186 DO jj = 1, jpj 187 DO ji = 1, jpi 188 ! 189 IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 190 191 ! find which category (jl0) the input ice thickness falls into 192 jl0 = jpl 193 DO jl = 1, jpl 194 IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 195 jl0 = jl 196 CYCLE 197 ENDIF 198 END DO 199 ! 200 itest(:) = 0 201 i_fill = jpl + 1 !------------------------------------ 202 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 203 ! !------------------------------------ 204 i_fill = i_fill - 1 205 ! 206 zh_i_ini(ji,jj,:) = 0._wp 207 za_i_ini(ji,jj,:) = 0._wp 208 itest(:) = 0 209 ! 210 IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1 211 zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj) 212 za_i_ini(ji,jj,1) = zat_i_ini(ji,jj) 213 ELSE !-- case ice is thicker: fill categories >1 214 ! thickness 215 DO jl = 1, i_fill-1 216 zh_i_ini(ji,jj,jl) = hi_mean(jl) 217 END DO 218 ! 219 ! concentration 220 za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 221 DO jl = 1, i_fill - 1 222 IF( jl /= jl0 )THEN 223 zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) ) 224 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 225 ENDIF 226 END DO 227 228 ! last category 229 za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) ) 230 zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) ) 231 zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 232 233 ! correction if concentration of upper cat is greater than lower cat 234 ! (it should be a gaussian around jl0 but sometimes it is not) 235 IF ( jl0 /= jpl ) THEN 236 DO jl = jpl, jl0+1, -1 237 IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN 238 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl) 239 zh_i_ini(ji,jj,jl ) = 0._wp 240 za_i_ini(ji,jj,jl ) = 0._wp 241 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1) & 242 & + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 ) 243 END IF 244 ENDDO 245 ENDIF 246 ! 247 ENDIF 248 ! 249 ! Compatibility tests 250 zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) ) ! Test 1: area conservation 251 IF ( zconv < epsi06 ) itest(1) = 1 252 ! 253 zconv = ABS( zat_i_ini(ji,jj) * zht_i_ini(ji,jj) & ! Test 2: volume conservation 254 & - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 255 IF ( zconv < epsi06 ) itest(2) = 1 256 ! 257 IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ? 258 ! 259 itest(4) = 1 260 DO jl = 1, i_fill 261 IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations 262 END DO 263 ! !---------------------------- 264 END DO ! end iteration on categories 265 ! !---------------------------- 266 IF( lwp .AND. SUM(itest) /= 4 ) THEN 267 WRITE(numout,*) 268 WRITE(numout,*) ' !!!! ALERT itest is not equal to 4 !!! ' 269 WRITE(numout,*) ' !!!! Something is wrong in the SI3 initialization procedure ' 270 WRITE(numout,*) 271 WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:) 272 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 273 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 274 ENDIF 275 ! 276 ENDIF 277 ! 278 END DO 279 END DO 280 ENDIF 281 282 !--------------------------------------------------------------------- 283 ! 4) Fill in sea ice arrays 284 !--------------------------------------------------------------------- 285 ! 286 ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 287 DO jl = 1, jpl ! loop over categories 288 DO jj = 1, jpj 289 DO ji = 1, jpi 290 a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini(ji,jj,jl) ! concentration 291 h_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(ji,jj,jl) ! ice thickness 292 s_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) ! salinity 293 o_i(ji,jj,jl) = 0._wp ! age (0 day) 294 t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 295 ! 296 IF( zht_i_ini(ji,jj) > 0._wp )THEN 297 h_s(ji,jj,jl)= h_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) ) ! snow depth 298 ELSE 299 h_s(ji,jj,jl)= 0._wp 300 ENDIF 301 ! 302 ! This case below should not be used if (h_s/h_i) is ok in namelist 303 ! In case snow load is in excess that would lead to transformation from snow to ice 304 ! Then, transfer the snow excess into the ice (different from icethd_dh) 305 zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 ) 306 ! recompute h_i, h_s avoiding out of bounds values 307 h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 308 h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi * r1_rhos ) 309 ! 310 ! ice volume, salt content, age content 311 v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume 312 v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume 313 sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 314 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content 315 END DO 316 END DO 133 ! ice and snw temperatures 134 DO jl = 1, jpl 135 DO jk = 1, nlay_i 136 t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) 317 137 END DO 318 !319 IF( nn_icesal /= 2 ) THEN ! for constant salinity in time320 CALL ice_var_salprof321 sv_i = s_i * v_i322 ENDIF323 !324 ! Snow temperature and heat content325 138 DO jk = 1, nlay_s 326 DO jl = 1, jpl ! loop over categories 327 DO jj = 1, jpj 328 DO ji = 1, jpi 329 t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 330 ! Snow energy of melting 331 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 332 ! 333 ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 334 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 335 END DO 336 END DO 337 END DO 139 t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) 338 140 END DO 339 !340 ! Ice salinity, temperature and heat content341 DO jk = 1, nlay_i342 DO jl = 1, jpl ! loop over categories343 DO jj = 1, jpj344 DO ji = 1, jpi345 t_i (ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0346 sz_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin347 ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K348 !349 ! heat content per unit volume350 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) &351 & + rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 ) ) &352 & - rcp * ( ztmelts - rt0 ) )353 !354 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2355 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i356 END DO357 END DO358 END DO359 END DO360 !361 tn_ice (:,:,:) = t_su (:,:,:)362 t1_ice (:,:,:) = t_i (:,:,1,:) ! initialisation of 1st layer temp for coupled simu363 364 ! Melt pond volume and fraction365 IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN ; zfac = 1._wp366 ELSE ; zfac = 0._wp367 ENDIF368 DO jl = 1, jpl369 a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac370 h_ip (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac371 END DO372 a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:)373 v_ip(:,:,:) = h_ip (:,:,:) * a_ip(:,:,:)374 !375 ELSE ! if ln_iceini=false376 a_i (:,:,:) = 0._wp377 v_i (:,:,:) = 0._wp378 v_s (:,:,:) = 0._wp379 sv_i (:,:,:) = 0._wp380 oa_i (:,:,:) = 0._wp381 h_i (:,:,:) = 0._wp382 h_s (:,:,:) = 0._wp383 s_i (:,:,:) = 0._wp384 o_i (:,:,:) = 0._wp385 !386 e_i(:,:,:,:) = 0._wp387 e_s(:,:,:,:) = 0._wp388 !389 DO jl = 1, jpl390 DO jk = 1, nlay_i391 t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)392 END DO393 DO jk = 1, nlay_s394 t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)395 END DO396 END DO397 398 tn_ice (:,:,:) = t_i (:,:,1,:)399 t1_ice (:,:,:) = t_i (:,:,1,:) ! initialisation of 1st layer temp for coupled simu400 401 a_ip(:,:,:) = 0._wp402 v_ip(:,:,:) = 0._wp403 a_ip_frac(:,:,:) = 0._wp404 h_ip (:,:,:) = 0._wp405 !406 ENDIF ! ln_iceini407 !408 at_i (:,:) = 0.0_wp409 DO jl = 1, jpl410 at_i (:,:) = at_i (:,:) + a_i (:,:,jl)411 141 END DO 412 142 ! 413 ! --- set ice velocities --- ! 143 ! specific temperatures for coupled runs 144 tn_ice (:,:,:) = t_i (:,:,1,:) 145 t1_ice (:,:,:) = t_i (:,:,1,:) 146 147 ! heat contents 148 e_i (:,:,:,:) = 0._wp 149 e_s (:,:,:,:) = 0._wp 150 151 ! general fields 152 a_i (:,:,:) = 0._wp 153 v_i (:,:,:) = 0._wp 154 v_s (:,:,:) = 0._wp 155 sv_i(:,:,:) = 0._wp 156 oa_i(:,:,:) = 0._wp 157 ! 158 h_i (:,:,:) = 0._wp 159 h_s (:,:,:) = 0._wp 160 s_i (:,:,:) = 0._wp 161 o_i (:,:,:) = 0._wp 162 ! 163 ! melt ponds 164 a_ip (:,:,:) = 0._wp 165 v_ip (:,:,:) = 0._wp 166 a_ip_frac(:,:,:) = 0._wp 167 h_ip (:,:,:) = 0._wp 168 ! 169 ! ice velocities 414 170 u_ice (:,:) = 0._wp 415 171 v_ice (:,:) = 0._wp 416 ! fields needed for ice_dyn_adv_umx 417 l_split_advumx(1) = .FALSE. 172 ! 173 !------------------------------------------------------------------------ 174 ! 2) overwrite some of the fields with namelist parameters or netcdf file 175 !------------------------------------------------------------------------ 176 177 178 IF( ln_iceini ) THEN 179 ! !---------------! 180 181 IF( Agrif_Root() ) THEN 182 183 IF( ln_iceini_file )THEN ! Read a file ! 184 ! !---------------! 185 WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp 186 ELSEWHERE ; zswitch(:,:) = 0._wp 187 END WHERE 188 ! 189 CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step 190 ! 191 ! -- mandatory fields -- ! 192 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1) 193 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1) 194 zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1) 195 196 ! -- optional fields -- ! 197 ! if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature) 198 ! 199 ! ice salinity 200 IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) & 201 & si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 202 ! 203 ! temperatures 204 IF ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. & 205 & TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN 206 si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 207 si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 208 si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 209 ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2 210 si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 ) 211 ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2 212 si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 ) 213 ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_su, set T_su = T_s 214 si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1) 215 ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_su, set T_su = T_i 216 si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 217 ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_s, set T_s = T_su 218 si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1) 219 ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_s, set T_s = T_i 220 si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 221 ENDIF 222 ! 223 ! pond concentration 224 IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 225 & si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc. 226 & * si(jp_ati)%fnow(:,:,1) 227 ! 228 ! pond depth 229 IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 230 & si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 231 ! 232 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1) 233 ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1) 234 zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1) 235 ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1) 236 zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1) 237 zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1) 238 ! 239 ! change the switch for the following 240 WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 241 ELSEWHERE ; zswitch(:,:) = 0._wp 242 END WHERE 243 244 ! !---------------! 245 ELSE ! Read namelist ! 246 ! !---------------! 247 ! no ice if (sst - Tfreez) >= thresold 248 WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 249 ELSEWHERE ; zswitch(:,:) = tmask(:,:,1) 250 END WHERE 251 ! 252 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 253 WHERE( ff_t(:,:) >= 0._wp ) 254 zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:) 255 zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 256 zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 257 zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 258 ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 259 zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 260 ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 261 zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 262 zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 263 ELSEWHERE 264 zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 265 zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 266 zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 267 zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 268 ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 269 zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:) 270 ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) 271 zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 272 zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 273 END WHERE 274 ! 275 ENDIF 276 277 278 279 ! make sure ponds = 0 if no ponds scheme 280 IF ( .NOT.ln_pnd ) THEN 281 zapnd_ini(:,:) = 0._wp 282 zhpnd_ini(:,:) = 0._wp 283 ENDIF 284 285 !-------------! 286 ! fill fields ! 287 !-------------! 288 ! select ice covered grid points 289 npti = 0 ; nptidx(:) = 0 290 DO_2D( 1, 1, 1, 1 ) 291 IF ( zht_i_ini(ji,jj) > 0._wp ) THEN 292 npti = npti + 1 293 nptidx(npti) = (jj - 1) * jpi + ji 294 ENDIF 295 END_2D 296 297 ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj) 298 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini ) 299 CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini ) 300 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini ) 301 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini ) 302 CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini ) 303 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini ) 304 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini ) 305 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini ) 306 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini ) 307 308 ! allocate temporary arrays 309 ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 310 & zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 311 312 ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 313 CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), & 314 & zhi_2d , zhs_2d , zai_2d , & 315 & t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), & 316 & zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d ) 317 318 ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 319 DO jl = 1, jpl 320 zti_3d(:,:,jl) = rt0 * tmask(:,:,1) 321 zts_3d(:,:,jl) = rt0 * tmask(:,:,1) 322 END DO 323 CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i ) 324 CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s ) 325 CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i ) 326 CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d ) 327 CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d ) 328 CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su ) 329 CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i ) 330 CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip ) 331 CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip ) 332 333 ! deallocate temporary arrays 334 DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 335 & zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d ) 336 337 ! calculate extensive and intensive variables 338 CALL ice_var_salprof ! for sz_i 339 DO jl = 1, jpl 340 DO_2D( 1, 1, 1, 1 ) 341 v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) 342 v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) 343 sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl) 344 END_2D 345 END DO 346 ! 347 DO jl = 1, jpl 348 DO_3D( 1, 1, 1, 1, 1, nlay_s ) 349 t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl) 350 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * & 351 & rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 352 END_3D 353 END DO 354 ! 355 DO jl = 1, jpl 356 DO_3D( 1, 1, 1, 1, 1, nlay_i ) 357 t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 358 ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K 359 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & 360 & rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + & 361 & rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) & 362 & - rcp * ( ztmelts - rt0 ) ) 363 END_3D 364 END DO 365 366 ! Melt ponds 367 WHERE( a_i > epsi10 ) 368 a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 369 ELSEWHERE 370 a_ip_frac(:,:,:) = 0._wp 371 END WHERE 372 v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 373 374 ! specific temperatures for coupled runs 375 tn_ice(:,:,:) = t_su(:,:,:) 376 t1_ice(:,:,:) = t_i (:,:,1,:) 377 ! 378 379 #if defined key_agrif 380 ELSE 381 382 Agrif_SpecialValue = -9999. 383 Agrif_UseSpecialValue = .TRUE. 384 CALL Agrif_init_variable(tra_iceini_id,procname=interp_tra_ice) 385 use_sign_north = .TRUE. 386 sign_north = -1. 387 CALL Agrif_init_variable(u_iceini_id ,procname=interp_u_ice) 388 CALL Agrif_init_variable(v_iceini_id ,procname=interp_v_ice) 389 Agrif_SpecialValue = 0._wp 390 use_sign_north = .FALSE. 391 Agrif_UseSpecialValue = .FALSE. 392 ! lbc ???? 393 ! Here we know : a_i, v_i, v_s, sv_i, oa_i, a_ip, v_ip, t_su, e_s, e_i 394 CALL ice_var_glo2eqv 395 CALL ice_var_zapsmall 396 CALL ice_var_agg(2) 397 398 ! Melt ponds 399 WHERE( a_i > epsi10 ) 400 a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 401 ELSEWHERE 402 a_ip_frac(:,:,:) = 0._wp 403 END WHERE 404 WHERE( a_ip > 0._wp ) ! ??????? 405 h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 406 ELSEWHERE 407 h_ip(:,:,:) = 0._wp 408 END WHERE 409 410 tn_ice(:,:,:) = t_su(:,:,:) 411 t1_ice(:,:,:) = t_i (:,:,1,:) 412 #endif 413 ENDIF ! Agrif_Root 414 ENDIF ! ln_iceini 415 ! 416 at_i(:,:) = SUM( a_i, dim=3 ) 418 417 ! 419 418 !---------------------------------------------- 420 ! 5) Snow-ice mass (case ice is fully embedded)419 ! 3) Snow-ice mass (case ice is fully embedded) 421 420 !---------------------------------------------- 422 421 snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass … … 425 424 IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area 426 425 ! 427 ssh n(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0428 ssh b(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0426 ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rho0 427 ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rho0 429 428 ! 430 IF( .NOT.ln_linssh ) THEN 431 ! 432 WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:) 433 ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE 434 ! 435 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 436 e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:) 437 e3t_b(:,:,jk) = e3t_n(:,:,jk) 438 e3t_a(:,:,jk) = e3t_n(:,:,jk) 439 END DO 440 ! 441 ! Reconstruction of all vertical scale factors at now and before time-steps 442 ! ========================================================================= 443 ! Horizontal scale factor interpolations 444 ! -------------------------------------- 445 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 446 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 447 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 448 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 449 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 450 ! Vertical scale factor interpolations 451 ! ------------------------------------ 452 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 453 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 454 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 455 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 456 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 457 ! t- and w- points depth 458 ! ---------------------- 459 !!gm not sure of that.... 460 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 461 gdepw_n(:,:,1) = 0.0_wp 462 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 463 DO jk = 2, jpk 464 gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk ) 465 gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) 466 gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:) 467 END DO 468 ENDIF 429 IF( .NOT.ln_linssh ) CALL dom_vvl_zgr( Kbb, Kmm, Kaa ) ! interpolation scale factor, depth and water column 430 ! !!st 431 ! IF( .NOT.ln_linssh ) THEN 432 ! ! 433 ! WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:) 434 ! ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE 435 ! ! 436 ! DO jk = 1,jpkm1 ! adjust initial vertical scale factors 437 ! e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:) 438 ! e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 439 ! e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm) 440 ! END DO 441 ! ! 442 ! ! Reconstruction of all vertical scale factors at now and before time-steps 443 ! ! ========================================================================= 444 ! ! Horizontal scale factor interpolations 445 ! ! -------------------------------------- 446 ! CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 447 ! CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 448 ! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 449 ! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 450 ! CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 451 ! ! Vertical scale factor interpolations 452 ! ! ------------------------------------ 453 ! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 454 ! CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 455 ! CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 456 ! CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 457 ! CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 458 ! ! t- and w- points depth 459 ! ! ---------------------- 460 ! !!gm not sure of that.... 461 ! gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 462 ! gdepw(:,:,1,Kmm) = 0.0_wp 463 ! gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 464 ! DO jk = 2, jpk 465 ! gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk ,Kmm) 466 ! gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm) 467 ! gde3w(:,:,jk) = gdept(:,:,jk ,Kmm) - ssh (:,:,Kmm) 468 ! END DO 469 ! ENDIF 469 470 ENDIF 470 471 471 472 !------------------------------------ 472 ! 6) store fields at before time-step473 ! 4) store fields at before time-step 473 474 !------------------------------------ 474 475 ! it is only necessary for the 1st interpolation by Agrif … … 487 488 !!clem: output of initial state should be written here but it is impossible because 488 489 !! the ocean and ice are in the same file 489 !! CALL dia_wri_state( 'output.init' )490 !! CALL dia_wri_state( Kmm, 'output.init' ) 490 491 ! 491 492 END SUBROUTINE ice_istate … … 504 505 !! 505 506 !!----------------------------------------------------------------------------- 506 INTEGER :: ios ! Local integer output status for namelist read507 INTEGER :: ifpr, ierror 507 INTEGER :: ios, ifpr, ierror ! Local integers 508 508 509 ! 509 510 CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files 510 TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_ tsu, sn_tmi, sn_smi511 TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd 511 512 TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read 512 513 ! 513 NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, & 514 & rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, & 515 & rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s, & 516 & sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir 514 NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, & 515 & rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, & 516 & rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, & 517 & rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, & 518 & rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, & 519 & sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir 517 520 !!----------------------------------------------------------------------------- 518 521 ! 519 REWIND( numnam_ice_ref ) ! Namelist namini in reference namelist : Ice initial state520 522 READ ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901) 521 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist', lwp ) 522 REWIND( numnam_ice_cfg ) ! Namelist namini in configuration namelist : Ice initial state 523 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' ) 523 524 READ ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 ) 524 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' , lwp)525 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' ) 525 526 IF(lwm) WRITE ( numoni, namini ) 526 527 ! 527 528 slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts 528 slf_i(jp_ati) = sn_ati ; slf_i(jp_tsu) = sn_tsu 529 slf_i(jp_tmi) = sn_tmi ; slf_i(jp_smi) = sn_smi 529 slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi 530 slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms 531 slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd 530 532 ! 531 533 IF(lwp) THEN ! control print … … 534 536 WRITE(numout,*) '~~~~~~~~~~~~~~~' 535 537 WRITE(numout,*) ' Namelist namini:' 536 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini 537 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file 538 WRITE(numout,*) ' max delta ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst 539 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n 540 WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s 541 WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n 542 WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s 543 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n 544 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s 545 WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n 546 WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s 547 WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n 548 WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s 538 WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini 539 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file 540 WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst 541 IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN 542 WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s 543 WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s 544 WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s 545 WRITE(numout,*) ' initial ice salinity in the north-south rn_smi_ini = ', rn_smi_ini_n,rn_smi_ini_s 546 WRITE(numout,*) ' initial surf temperat in the north-south rn_tsu_ini = ', rn_tsu_ini_n,rn_tsu_ini_s 547 WRITE(numout,*) ' initial ice temperat in the north-south rn_tmi_ini = ', rn_tmi_ini_n,rn_tmi_ini_s 548 WRITE(numout,*) ' initial snw temperat in the north-south rn_tms_ini = ', rn_tms_ini_n,rn_tms_ini_s 549 WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s 550 WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s 551 ENDIF 549 552 ENDIF 550 553 ! … … 554 557 ALLOCATE( si(jpfldi), STAT=ierror ) 555 558 IF( ierror > 0 ) THEN 556 CALL ctl_stop( ' Ice_ini in iceistate: unable to allocate si structure' ) ; RETURN559 CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' ) ; RETURN 557 560 ENDIF 558 561 ! 559 562 DO ifpr = 1, jpfldi 560 563 ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) 561 ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )564 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 562 565 END DO 563 566 ! 564 567 ! fill si with slf_i and control print 565 CALL fld_fill( si, slf_i, cn_dir, 'ice_istate', 'ice istate ini', 'numnam_ice' ) 566 ! 567 CALL fld_read( nit000, 1, si ) ! input fields provided at the current time-step 568 CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' ) 568 569 ! 569 570 ENDIF 570 571 ! 572 IF( .NOT.ln_pnd ) THEN 573 rn_apd_ini_n = 0. ; rn_apd_ini_s = 0. 574 rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0. 575 CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' ) 576 ENDIF 577 ! 571 578 END SUBROUTINE ice_istate_init 572 579
Note: See TracChangeset
for help on using the changeset viewer.