Changeset 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/iceistate.F90
- Timestamp:
- 2019-10-29T11:41:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/iceistate.F90
r10998 r11822 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 … … 36 38 PUBLIC ice_istate ! called by icestp.F90 37 39 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 40 ! 48 41 ! !! ** 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 42 LOGICAL, PUBLIC :: ln_iceini !: Ice initialization or not 43 LOGICAL, PUBLIC :: ln_iceini_file !: Ice initialization from 2D netcdf file 44 REAL(wp) :: rn_thres_sst 45 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 46 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 47 REAL(wp) :: rn_apd_ini_n, rn_hpd_ini_n 48 REAL(wp) :: rn_apd_ini_s, rn_hpd_ini_s 49 ! 50 ! ! if ln_iceini_file = T 51 INTEGER , PARAMETER :: jpfldi = 9 ! maximum number of files to read 52 INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m) 53 INTEGER , PARAMETER :: jp_hts = 2 ! index of snw thickness (m) 54 INTEGER , PARAMETER :: jp_ati = 3 ! index of ice fraction (-) 55 INTEGER , PARAMETER :: jp_smi = 4 ! index of ice salinity (g/kg) 56 INTEGER , PARAMETER :: jp_tmi = 5 ! index of ice temperature (K) 57 INTEGER , PARAMETER :: jp_tsu = 6 ! index of ice surface temp (K) 58 INTEGER , PARAMETER :: jp_tms = 7 ! index of snw temperature (K) 59 INTEGER , PARAMETER :: jp_apd = 8 ! index of pnd fraction (-) 60 INTEGER , PARAMETER :: jp_hpd = 9 ! index of pnd depth (m) 61 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read) 62 ! 63 63 !!---------------------------------------------------------------------- 64 64 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 68 68 CONTAINS 69 69 70 SUBROUTINE ice_istate( Kbb, Kmm, Kaa )70 SUBROUTINE ice_istate( kt, Kbb, Kmm, Kaa ) 71 71 !!------------------------------------------------------------------- 72 72 !! *** ROUTINE ice_istate *** … … 87 87 !! 88 88 !! ** 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?)89 !! where there is no ice 90 90 !!-------------------------------------------------------------------- 91 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 91 INTEGER, INTENT(in) :: kt ! time step 92 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 92 93 ! 93 94 INTEGER :: ji, jj, jk, jl ! dummy loop indices 94 INTEGER :: i_hemis, i_fill, jl0 ! local integers 95 REAL(wp) :: ztmelts, zdh 96 REAL(wp) :: zarg, zV, zconv, zdv, zfac 95 REAL(wp) :: ztmelts 97 96 INTEGER , DIMENSION(4) :: itest 98 97 REAL(wp), DIMENSION(jpi,jpj) :: z2d 99 98 REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator 100 REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file 101 REAL(wp), DIMENSION(jpi,jpj) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 102 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini , za_i_ini !data by cattegories to fill 99 REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file 100 REAL(wp), DIMENSION(jpi,jpj) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 101 REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini !data from namelist or nc file 102 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !temporary arrays 103 !! 104 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d 103 105 !-------------------------------------------------------------------- 104 106 … … 107 109 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 108 110 109 !-------------------------------------------------------------------- 110 ! 1) Set surface and bottom temperatures to initial values 111 !-------------------------------------------------------------------- 112 ! 113 ! init surface temperature 111 !--------------------------- 112 ! 1) 1st init. of the fields 113 !--------------------------- 114 ! 115 ! basal temperature (considered at freezing point) [Kelvin] 116 CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 117 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 118 ! 119 ! surface temperature and conductivity 114 120 DO jl = 1, jpl 115 121 t_su (:,:,jl) = rt0 * tmask(:,:,1) ! temp at the surface … … 117 123 END DO 118 124 ! 119 ! init basal temperature (considered at freezing point) [Kelvin] 120 CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 121 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 122 125 ! ice and snw temperatures 126 DO jl = 1, jpl 127 DO jk = 1, nlay_i 128 t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) 129 END DO 130 DO jk = 1, nlay_s 131 t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) 132 END DO 133 END DO 134 ! 135 ! specific temperatures for coupled runs 136 tn_ice (:,:,:) = t_i (:,:,1,:) 137 t1_ice (:,:,:) = t_i (:,:,1,:) 138 139 ! heat contents 140 e_i (:,:,:,:) = 0._wp 141 e_s (:,:,:,:) = 0._wp 142 143 ! general fields 144 a_i (:,:,:) = 0._wp 145 v_i (:,:,:) = 0._wp 146 v_s (:,:,:) = 0._wp 147 sv_i(:,:,:) = 0._wp 148 oa_i(:,:,:) = 0._wp 149 ! 150 h_i (:,:,:) = 0._wp 151 h_s (:,:,:) = 0._wp 152 s_i (:,:,:) = 0._wp 153 o_i (:,:,:) = 0._wp 154 ! 155 ! melt ponds 156 a_ip (:,:,:) = 0._wp 157 v_ip (:,:,:) = 0._wp 158 a_ip_frac(:,:,:) = 0._wp 159 h_ip (:,:,:) = 0._wp 160 ! 161 ! ice velocities 162 u_ice (:,:) = 0._wp 163 v_ice (:,:) = 0._wp 164 ! 165 !------------------------------------------------------------------------ 166 ! 2) overwrite some of the fields with namelist parameters or netcdf file 167 !------------------------------------------------------------------------ 123 168 IF( ln_iceini ) THEN 124 !-----------------------------------------------------------125 ! 2) Compute or read sea ice variables ===> single category126 !-----------------------------------------------------------127 !128 169 ! !---------------! 129 170 IF( ln_iceini_file )THEN ! Read a file ! 130 171 ! !---------------! 131 ! 132 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) 133 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) 134 zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) 135 zts_u_ini(:,:) = si(jp_tsu)%fnow(:,:,1) 136 ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) 137 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) 138 ! 139 WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 140 ELSEWHERE ; zswitch(:,:) = 0._wp 172 WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp 173 ELSEWHERE ; zswitch(:,:) = 0._wp 141 174 END WHERE 142 zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) 143 ! 175 ! 176 CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step 177 ! 178 ! -- mandatory fields -- ! 179 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) 180 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) 181 zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) 182 183 ! -- optional fields -- ! 184 ! if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature) 185 ! 186 ! ice salinity 187 IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) & 188 & si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 189 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) 190 ! 191 ! ice temperature 192 IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) & 193 & si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 194 ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) 195 ! 196 ! surface temperature => set to ice temperature if it exists 197 IF ( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN 198 si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 199 ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN 200 si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 201 ENDIF 202 zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) 203 ! 204 ! snow temperature => set to ice temperature if it exists 205 IF ( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN 206 si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 207 ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN 208 si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 209 ENDIF 210 ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) 211 ! 212 ! pond concentration 213 IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 214 & 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. 215 & * si(jp_ati)%fnow(:,:,1) 216 zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) 217 ! 218 ! pond depth 219 IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 220 & si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 221 zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) 222 ! 223 ! change the switch for the following 224 WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 225 ELSEWHERE ; zswitch(:,:) = 0._wp 226 END WHERE 144 227 ! !---------------! 145 228 ELSE ! Read namelist ! 146 229 ! !---------------! 147 ! no ice if sst <= t-freez + ttest230 ! no ice if (sst - Tfreez) >= thresold 148 231 WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 149 232 ELSEWHERE ; zswitch(:,:) = tmask(:,:,1) … … 155 238 zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 156 239 zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 157 zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)158 240 zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 159 241 ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 242 zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 243 ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 244 zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 245 zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 160 246 ELSEWHERE 161 247 zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 162 248 zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 163 249 zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 164 zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)165 250 zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 166 251 ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 252 zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:) 253 ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) 254 zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 255 zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 167 256 END WHERE 168 zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) 169 ! 257 ! 258 ENDIF 259 260 ! make sure ponds = 0 if no ponds scheme 261 IF ( .NOT.ln_pnd ) THEN 262 zapnd_ini(:,:) = 0._wp 263 zhpnd_ini(:,:) = 0._wp 170 264 ENDIF 171 265 172 !------------------------------------------------------------------ 173 ! 3) Distribute ice concentration and thickness into the categories 174 !------------------------------------------------------------------ 175 ! a gaussian distribution for ice concentration is used 176 ! then we check whether the distribution fullfills 177 ! volume and area conservation, positivity and ice categories bounds 178 179 IF( jpl == 1 ) THEN 180 ! 181 zh_i_ini(:,:,1) = zht_i_ini(:,:) 182 za_i_ini(:,:,1) = zat_i_ini(:,:) 183 ! 184 ELSE 185 zh_i_ini(:,:,:) = 0._wp 186 za_i_ini(:,:,:) = 0._wp 187 ! 266 !-------------! 267 ! fill fields ! 268 !-------------! 269 ! select ice covered grid points 270 npti = 0 ; nptidx(:) = 0 271 DO jj = 1, jpj 272 DO ji = 1, jpi 273 IF ( zht_i_ini(ji,jj) > 0._wp ) THEN 274 npti = npti + 1 275 nptidx(npti) = (jj - 1) * jpi + ji 276 ENDIF 277 END DO 278 END DO 279 280 ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj) 281 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini ) 282 CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini ) 283 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini ) 284 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini ) 285 CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini ) 286 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini ) 287 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini ) 288 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini ) 289 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini ) 290 291 ! allocate temporary arrays 292 ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 293 & zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 294 295 ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 296 CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), & 297 & zhi_2d , zhs_2d , zai_2d , & 298 & 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), & 299 & zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d ) 300 301 ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 302 DO jl = 1, jpl 303 zti_3d(:,:,jl) = rt0 * tmask(:,:,1) 304 zts_3d(:,:,jl) = rt0 * tmask(:,:,1) 305 END DO 306 CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i ) 307 CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s ) 308 CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i ) 309 CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d ) 310 CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d ) 311 CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su ) 312 CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i ) 313 CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip ) 314 CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip ) 315 316 ! deallocate temporary arrays 317 DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 318 & zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d ) 319 320 ! calculate extensive and intensive variables 321 CALL ice_var_salprof ! for sz_i 322 DO jl = 1, jpl 188 323 DO jj = 1, jpj 189 324 DO ji = 1, jpi 190 ! 191 IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 192 193 ! find which category (jl0) the input ice thickness falls into 194 jl0 = jpl 195 DO jl = 1, jpl 196 IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 197 jl0 = jl 198 CYCLE 199 ENDIF 200 END DO 201 ! 202 itest(:) = 0 203 i_fill = jpl + 1 !------------------------------------ 204 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 205 ! !------------------------------------ 206 i_fill = i_fill - 1 207 ! 208 zh_i_ini(ji,jj,:) = 0._wp 209 za_i_ini(ji,jj,:) = 0._wp 210 itest(:) = 0 211 ! 212 IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1 213 zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj) 214 za_i_ini(ji,jj,1) = zat_i_ini(ji,jj) 215 ELSE !-- case ice is thicker: fill categories >1 216 ! thickness 217 DO jl = 1, i_fill-1 218 zh_i_ini(ji,jj,jl) = hi_mean(jl) 219 END DO 220 ! 221 ! concentration 222 za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 223 DO jl = 1, i_fill - 1 224 IF( jl /= jl0 )THEN 225 zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) ) 226 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 227 ENDIF 228 END DO 229 230 ! last category 231 za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) ) 232 zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) ) 233 zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 234 235 ! correction if concentration of upper cat is greater than lower cat 236 ! (it should be a gaussian around jl0 but sometimes it is not) 237 IF ( jl0 /= jpl ) THEN 238 DO jl = jpl, jl0+1, -1 239 IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN 240 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl) 241 zh_i_ini(ji,jj,jl ) = 0._wp 242 za_i_ini(ji,jj,jl ) = 0._wp 243 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1) & 244 & + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 ) 245 END IF 246 ENDDO 247 ENDIF 248 ! 249 ENDIF 250 ! 251 ! Compatibility tests 252 zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) ) ! Test 1: area conservation 253 IF ( zconv < epsi06 ) itest(1) = 1 254 ! 255 zconv = ABS( zat_i_ini(ji,jj) * zht_i_ini(ji,jj) & ! Test 2: volume conservation 256 & - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 257 IF ( zconv < epsi06 ) itest(2) = 1 258 ! 259 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 ? 260 ! 261 itest(4) = 1 262 DO jl = 1, i_fill 263 IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations 264 END DO 265 ! !---------------------------- 266 END DO ! end iteration on categories 267 ! !---------------------------- 268 IF( lwp .AND. SUM(itest) /= 4 ) THEN 269 WRITE(numout,*) 270 WRITE(numout,*) ' !!!! ALERT itest is not equal to 4 !!! ' 271 WRITE(numout,*) ' !!!! Something is wrong in the SI3 initialization procedure ' 272 WRITE(numout,*) 273 WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:) 274 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 275 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 276 ENDIF 277 ! 278 ENDIF 279 ! 325 v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) 326 v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) 327 sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl) 280 328 END DO 281 329 END DO 282 ENDIF 283 284 !--------------------------------------------------------------------- 285 ! 4) Fill in sea ice arrays 286 !--------------------------------------------------------------------- 287 ! 288 ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 289 DO jl = 1, jpl ! loop over categories 290 DO jj = 1, jpj 291 DO ji = 1, jpi 292 a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini(ji,jj,jl) ! concentration 293 h_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(ji,jj,jl) ! ice thickness 294 s_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) ! salinity 295 o_i(ji,jj,jl) = 0._wp ! age (0 day) 296 t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 297 ! 298 IF( zht_i_ini(ji,jj) > 0._wp )THEN 299 h_s(ji,jj,jl)= h_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) ) ! snow depth 300 ELSE 301 h_s(ji,jj,jl)= 0._wp 302 ENDIF 303 ! 304 ! This case below should not be used if (h_s/h_i) is ok in namelist 305 ! In case snow load is in excess that would lead to transformation from snow to ice 306 ! Then, transfer the snow excess into the ice (different from icethd_dh) 307 zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 ) 308 ! recompute h_i, h_s avoiding out of bounds values 309 h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 310 h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi * r1_rhos ) 311 ! 312 ! ice volume, salt content, age content 313 v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume 314 v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume 315 sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 316 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content 317 END DO 318 END DO 319 END DO 320 ! 321 IF( nn_icesal /= 2 ) THEN ! for constant salinity in time 322 CALL ice_var_salprof 323 sv_i = s_i * v_i 324 ENDIF 325 ! 326 ! Snow temperature and heat content 327 DO jk = 1, nlay_s 328 DO jl = 1, jpl ! loop over categories 330 END DO 331 ! 332 DO jl = 1, jpl 333 DO jk = 1, nlay_s 329 334 DO jj = 1, jpj 330 335 DO ji = 1, jpi 331 t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 332 ! Snow energy of melting 333 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 334 ! 335 ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 336 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 336 t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl) 337 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * & 338 & rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 337 339 END DO 338 340 END DO … … 340 342 END DO 341 343 ! 342 ! Ice salinity, temperature and heat content 343 DO jk = 1, nlay_i 344 DO jl = 1, jpl ! loop over categories 344 DO jl = 1, jpl 345 DO jk = 1, nlay_i 345 346 DO jj = 1, jpj 346 347 DO ji = 1, jpi 347 t_i (ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 348 sz_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin 349 ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 350 ! 351 ! heat content per unit volume 352 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) & 353 & + rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 ) ) & 354 & - rcp * ( ztmelts - rt0 ) ) 355 ! 356 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 357 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 348 t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 349 ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K 350 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & 351 & rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + & 352 & rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) & 353 & - rcp * ( ztmelts - rt0 ) ) 358 354 END DO 359 355 END DO 360 356 END DO 361 357 END DO 362 ! 363 tn_ice (:,:,:) = t_su (:,:,:) 364 t1_ice (:,:,:) = t_i (:,:,1,:) ! initialisation of 1st layer temp for coupled simu 365 366 ! Melt pond volume and fraction 367 IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN ; zfac = 1._wp 368 ELSE ; zfac = 0._wp 369 ENDIF 370 DO jl = 1, jpl 371 a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac 372 h_ip (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac 373 END DO 374 a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:) 375 v_ip(:,:,:) = h_ip (:,:,:) * a_ip(:,:,:) 376 ! 377 ELSE ! if ln_iceini=false 378 a_i (:,:,:) = 0._wp 379 v_i (:,:,:) = 0._wp 380 v_s (:,:,:) = 0._wp 381 sv_i (:,:,:) = 0._wp 382 oa_i (:,:,:) = 0._wp 383 h_i (:,:,:) = 0._wp 384 h_s (:,:,:) = 0._wp 385 s_i (:,:,:) = 0._wp 386 o_i (:,:,:) = 0._wp 387 ! 388 e_i(:,:,:,:) = 0._wp 389 e_s(:,:,:,:) = 0._wp 390 ! 391 DO jl = 1, jpl 392 DO jk = 1, nlay_i 393 t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) 394 END DO 395 DO jk = 1, nlay_s 396 t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) 397 END DO 398 END DO 399 400 tn_ice (:,:,:) = t_i (:,:,1,:) 401 t1_ice (:,:,:) = t_i (:,:,1,:) ! initialisation of 1st layer temp for coupled simu 402 403 a_ip(:,:,:) = 0._wp 404 v_ip(:,:,:) = 0._wp 405 a_ip_frac(:,:,:) = 0._wp 406 h_ip (:,:,:) = 0._wp 358 359 ! Melt ponds 360 WHERE( a_i > epsi10 ) 361 a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 362 ELSEWHERE 363 a_ip_frac(:,:,:) = 0._wp 364 END WHERE 365 v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 366 367 ! specific temperatures for coupled runs 368 tn_ice(:,:,:) = t_su(:,:,:) 369 t1_ice(:,:,:) = t_i (:,:,1,:) 407 370 ! 408 371 ENDIF ! ln_iceini 409 372 ! 410 at_i (:,:) = 0.0_wp 411 DO jl = 1, jpl 412 at_i (:,:) = at_i (:,:) + a_i (:,:,jl) 413 END DO 414 ! 415 ! --- set ice velocities --- ! 416 u_ice (:,:) = 0._wp 417 v_ice (:,:) = 0._wp 418 ! fields needed for ice_dyn_adv_umx 419 l_split_advumx(1) = .FALSE. 373 at_i(:,:) = SUM( a_i, dim=3 ) 420 374 ! 421 375 !---------------------------------------------- 422 ! 5) Snow-ice mass (case ice is fully embedded)376 ! 3) Snow-ice mass (case ice is fully embedded) 423 377 !---------------------------------------------- 424 378 snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass … … 472 426 473 427 !------------------------------------ 474 ! 6) store fields at before time-step428 ! 4) store fields at before time-step 475 429 !------------------------------------ 476 430 ! it is only necessary for the 1st interpolation by Agrif … … 506 460 !! 507 461 !!----------------------------------------------------------------------------- 508 INTEGER :: ji, jj 509 INTEGER :: ios, ierr, inum_ice ! Local integer output status for namelist read 462 INTEGER :: ios ! Local integer output status for namelist read 510 463 INTEGER :: ifpr, ierror 511 464 ! 512 465 CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files 513 TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_ tsu, sn_tmi, sn_smi466 TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd 514 467 TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read 515 468 ! 516 NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, & 517 & rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, & 518 & rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s, & 519 & sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir 469 NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, & 470 & rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, & 471 & rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, & 472 & rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, & 473 & rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, & 474 & sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir 520 475 !!----------------------------------------------------------------------------- 521 476 ! 522 477 REWIND( numnam_ice_ref ) ! Namelist namini in reference namelist : Ice initial state 523 478 READ ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901) 524 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' , lwp)479 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' ) 525 480 REWIND( numnam_ice_cfg ) ! Namelist namini in configuration namelist : Ice initial state 526 481 READ ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 ) 527 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' , lwp)482 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' ) 528 483 IF(lwm) WRITE ( numoni, namini ) 529 484 ! 530 485 slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts 531 slf_i(jp_ati) = sn_ati ; slf_i(jp_tsu) = sn_tsu 532 slf_i(jp_tmi) = sn_tmi ; slf_i(jp_smi) = sn_smi 486 slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi 487 slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms 488 slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd 533 489 ! 534 490 IF(lwp) THEN ! control print … … 537 493 WRITE(numout,*) '~~~~~~~~~~~~~~~' 538 494 WRITE(numout,*) ' Namelist namini:' 539 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini 540 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file 541 WRITE(numout,*) ' max delta ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst 542 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n 543 WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s 544 WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n 545 WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s 546 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n 547 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s 548 WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n 549 WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s 550 WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n 551 WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s 495 WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini 496 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file 497 WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst 498 IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN 499 WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s 500 WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s 501 WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s 502 WRITE(numout,*) ' initial ice salinity in the north-south rn_smi_ini = ', rn_smi_ini_n,rn_smi_ini_s 503 WRITE(numout,*) ' initial surf temperat in the north-south rn_tsu_ini = ', rn_tsu_ini_n,rn_tsu_ini_s 504 WRITE(numout,*) ' initial ice temperat in the north-south rn_tmi_ini = ', rn_tmi_ini_n,rn_tmi_ini_s 505 WRITE(numout,*) ' initial snw temperat in the north-south rn_tms_ini = ', rn_tms_ini_n,rn_tms_ini_s 506 WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s 507 WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s 508 ENDIF 552 509 ENDIF 553 510 ! … … 557 514 ALLOCATE( si(jpfldi), STAT=ierror ) 558 515 IF( ierror > 0 ) THEN 559 CALL ctl_stop( ' Ice_ini in iceistate: unable to allocate si structure' ) ; RETURN516 CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' ) ; RETURN 560 517 ENDIF 561 518 ! 562 519 DO ifpr = 1, jpfldi 563 520 ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) 564 ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )521 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 565 522 END DO 566 523 ! 567 524 ! fill si with slf_i and control print 568 CALL fld_fill( si, slf_i, cn_dir, 'ice_istate', 'ice istate ini', 'numnam_ice' ) 569 ! 570 CALL fld_read( nit000, 1, si ) ! input fields provided at the current time-step 571 ! 525 CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' ) 526 ! 527 ENDIF 528 ! 529 IF( .NOT.ln_pnd ) THEN 530 rn_apd_ini_n = 0. ; rn_apd_ini_s = 0. 531 rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0. 532 CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' ) 572 533 ENDIF 573 534 !
Note: See TracChangeset
for help on using the changeset viewer.
