Changeset 11573 for NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/ICE/iceistate.F90
- Timestamp:
- 2019-09-19T11:18:03+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/ICE/iceistate.F90
r11229 r11573 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 70 SUBROUTINE ice_istate( kt ) 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) :: kt ! time step 92 !! 91 93 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 94 REAL(wp) :: ztmelts 95 95 INTEGER , DIMENSION(4) :: itest 96 96 REAL(wp), DIMENSION(jpi,jpj) :: z2d 97 97 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 98 REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file 99 REAL(wp), DIMENSION(jpi,jpj) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 100 REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini !data from namelist or nc file 101 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !temporary arrays 102 !! 103 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d 101 104 !-------------------------------------------------------------------- 102 105 … … 105 108 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 106 109 107 !-------------------------------------------------------------------- 108 ! 1) Set surface and bottom temperatures to initial values 109 !-------------------------------------------------------------------- 110 ! 111 ! init surface temperature 110 !--------------------------- 111 ! 1) 1st init. of the fields 112 !--------------------------- 113 ! 114 ! basal temperature (considered at freezing point) [Kelvin] 115 CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 116 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 117 ! 118 ! surface temperature and conductivity 112 119 DO jl = 1, jpl 113 120 t_su (:,:,jl) = rt0 * tmask(:,:,1) ! temp at the surface … … 115 122 END DO 116 123 ! 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 124 ! ice and snw temperatures 125 DO jl = 1, jpl 126 DO jk = 1, nlay_i 127 t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) 128 END DO 129 DO jk = 1, nlay_s 130 t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) 131 END DO 132 END DO 133 ! 134 ! specific temperatures for coupled runs 135 tn_ice (:,:,:) = t_i (:,:,1,:) 136 t1_ice (:,:,:) = t_i (:,:,1,:) 137 138 ! heat contents 139 e_i (:,:,:,:) = 0._wp 140 e_s (:,:,:,:) = 0._wp 141 142 ! general fields 143 a_i (:,:,:) = 0._wp 144 v_i (:,:,:) = 0._wp 145 v_s (:,:,:) = 0._wp 146 sv_i(:,:,:) = 0._wp 147 oa_i(:,:,:) = 0._wp 148 ! 149 h_i (:,:,:) = 0._wp 150 h_s (:,:,:) = 0._wp 151 s_i (:,:,:) = 0._wp 152 o_i (:,:,:) = 0._wp 153 ! 154 ! melt ponds 155 a_ip (:,:,:) = 0._wp 156 v_ip (:,:,:) = 0._wp 157 a_ip_frac(:,:,:) = 0._wp 158 h_ip (:,:,:) = 0._wp 159 ! 160 ! ice velocities 161 u_ice (:,:) = 0._wp 162 v_ice (:,:) = 0._wp 163 ! 164 !------------------------------------------------------------------------ 165 ! 2) overwrite some of the fields with namelist parameters or netcdf file 166 !------------------------------------------------------------------------ 121 167 IF( ln_iceini ) THEN 122 !-----------------------------------------------------------123 ! 2) Compute or read sea ice variables ===> single category124 !-----------------------------------------------------------125 !126 168 ! !---------------! 127 169 IF( ln_iceini_file )THEN ! Read a file ! 128 170 ! !---------------! 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 171 WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp 172 ELSEWHERE ; zswitch(:,:) = 0._wp 139 173 END WHERE 140 zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) 141 ! 174 ! 175 CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step 176 ! 177 ! -- mandatory fields -- ! 178 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) 179 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) 180 zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) 181 182 ! -- optional fields -- ! 183 ! if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature) 184 ! 185 ! ice salinity 186 IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) & 187 & si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 188 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) 189 ! 190 ! ice temperature 191 IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) & 192 & si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 193 ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) 194 ! 195 ! surface temperature => set to ice temperature if it exists 196 IF ( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN 197 si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 198 ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN 199 si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 200 ENDIF 201 zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) 202 ! 203 ! snow temperature => set to ice temperature if it exists 204 IF ( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN 205 si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 206 ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN 207 si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 208 ENDIF 209 ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) 210 ! 211 ! pond concentration 212 IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 213 & 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. 214 & * si(jp_ati)%fnow(:,:,1) 215 zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) 216 ! 217 ! pond depth 218 IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 219 & si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 220 zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) 221 ! 222 ! change the switch for the following 223 WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 224 ELSEWHERE ; zswitch(:,:) = 0._wp 225 END WHERE 142 226 ! !---------------! 143 227 ELSE ! Read namelist ! 144 228 ! !---------------! 145 ! no ice if sst <= t-freez + ttest229 ! no ice if (sst - Tfreez) >= thresold 146 230 WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 147 231 ELSEWHERE ; zswitch(:,:) = tmask(:,:,1) … … 153 237 zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 154 238 zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 155 zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)156 239 zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 157 240 ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 241 zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 242 ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 243 zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 244 zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 158 245 ELSEWHERE 159 246 zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 160 247 zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 161 248 zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 162 zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)163 249 zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 164 250 ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 251 zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:) 252 ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) 253 zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 254 zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 165 255 END WHERE 166 zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) 167 ! 256 ! 257 ENDIF 258 259 ! make sure ponds = 0 if no ponds scheme 260 IF ( .NOT.ln_pnd ) THEN 261 zapnd_ini(:,:) = 0._wp 262 zhpnd_ini(:,:) = 0._wp 168 263 ENDIF 169 264 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 ! 265 !-------------! 266 ! fill fields ! 267 !-------------! 268 ! select ice covered grid points 269 npti = 0 ; nptidx(:) = 0 270 DO jj = 1, jpj 271 DO ji = 1, jpi 272 IF ( zht_i_ini(ji,jj) > 0._wp ) THEN 273 npti = npti + 1 274 nptidx(npti) = (jj - 1) * jpi + ji 275 ENDIF 276 END DO 277 END DO 278 279 ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj) 280 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini ) 281 CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini ) 282 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini ) 283 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini ) 284 CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini ) 285 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini ) 286 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini ) 287 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini ) 288 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini ) 289 290 ! allocate temporary arrays 291 ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 292 & zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 293 294 ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 295 CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), & 296 & zhi_2d , zhs_2d , zai_2d , & 297 & 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), & 298 & zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d ) 299 300 ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 301 DO jl = 1, jpl 302 zti_3d(:,:,jl) = rt0 * tmask(:,:,1) 303 zts_3d(:,:,jl) = rt0 * tmask(:,:,1) 304 END DO 305 CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i ) 306 CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s ) 307 CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i ) 308 CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d ) 309 CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d ) 310 CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su ) 311 CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i ) 312 CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip ) 313 CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip ) 314 315 ! deallocate temporary arrays 316 DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 317 & zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d ) 318 319 ! calculate extensive and intensive variables 320 CALL ice_var_salprof ! for sz_i 321 DO jl = 1, jpl 186 322 DO jj = 1, jpj 187 323 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 ! 324 v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) 325 v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) 326 sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl) 278 327 END DO 279 328 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 317 END DO 318 ! 319 IF( nn_icesal /= 2 ) THEN ! for constant salinity in time 320 CALL ice_var_salprof 321 sv_i = s_i * v_i 322 ENDIF 323 ! 324 ! Snow temperature and heat content 325 DO jk = 1, nlay_s 326 DO jl = 1, jpl ! loop over categories 329 END DO 330 ! 331 DO jl = 1, jpl 332 DO jk = 1, nlay_s 327 333 DO jj = 1, jpj 328 334 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 t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl) 336 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * & 337 & rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 335 338 END DO 336 339 END DO … … 338 341 END DO 339 342 ! 340 ! Ice salinity, temperature and heat content 341 DO jk = 1, nlay_i 342 DO jl = 1, jpl ! loop over categories 343 DO jl = 1, jpl 344 DO jk = 1, nlay_i 343 345 DO jj = 1, jpj 344 346 DO ji = 1, jpi 345 t_i (ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 346 sz_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin 347 ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 348 ! 349 ! heat content per unit volume 350 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/m2 355 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 347 t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 348 ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K 349 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & 350 & 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 ) ) 356 353 END DO 357 354 END DO 358 355 END DO 359 356 END DO 360 ! 361 tn_ice (:,:,:) = t_su (:,:,:) 362 t1_ice (:,:,:) = t_i (:,:,1,:) ! initialisation of 1st layer temp for coupled simu 363 364 ! Melt pond volume and fraction 365 IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN ; zfac = 1._wp 366 ELSE ; zfac = 0._wp 367 ENDIF 368 DO jl = 1, jpl 369 a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac 370 h_ip (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac 371 END DO 372 a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:) 373 v_ip(:,:,:) = h_ip (:,:,:) * a_ip(:,:,:) 374 ! 375 ELSE ! if ln_iceini=false 376 a_i (:,:,:) = 0._wp 377 v_i (:,:,:) = 0._wp 378 v_s (:,:,:) = 0._wp 379 sv_i (:,:,:) = 0._wp 380 oa_i (:,:,:) = 0._wp 381 h_i (:,:,:) = 0._wp 382 h_s (:,:,:) = 0._wp 383 s_i (:,:,:) = 0._wp 384 o_i (:,:,:) = 0._wp 385 ! 386 e_i(:,:,:,:) = 0._wp 387 e_s(:,:,:,:) = 0._wp 388 ! 389 DO jl = 1, jpl 390 DO jk = 1, nlay_i 391 t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) 392 END DO 393 DO jk = 1, nlay_s 394 t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) 395 END DO 396 END DO 397 398 tn_ice (:,:,:) = t_i (:,:,1,:) 399 t1_ice (:,:,:) = t_i (:,:,1,:) ! initialisation of 1st layer temp for coupled simu 400 401 a_ip(:,:,:) = 0._wp 402 v_ip(:,:,:) = 0._wp 403 a_ip_frac(:,:,:) = 0._wp 404 h_ip (:,:,:) = 0._wp 357 358 ! Melt ponds 359 WHERE( a_i > epsi10 ) 360 a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 361 ELSEWHERE 362 a_ip_frac(:,:,:) = 0._wp 363 END WHERE 364 v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 365 366 ! specific temperatures for coupled runs 367 tn_ice(:,:,:) = t_su(:,:,:) 368 t1_ice(:,:,:) = t_i (:,:,1,:) 405 369 ! 406 370 ENDIF ! ln_iceini 407 371 ! 408 at_i (:,:) = 0.0_wp 409 DO jl = 1, jpl 410 at_i (:,:) = at_i (:,:) + a_i (:,:,jl) 411 END DO 412 ! 413 ! --- set ice velocities --- ! 414 u_ice (:,:) = 0._wp 415 v_ice (:,:) = 0._wp 416 ! fields needed for ice_dyn_adv_umx 417 l_split_advumx(1) = .FALSE. 372 at_i(:,:) = SUM( a_i, dim=3 ) 418 373 ! 419 374 !---------------------------------------------- 420 ! 5) Snow-ice mass (case ice is fully embedded)375 ! 3) Snow-ice mass (case ice is fully embedded) 421 376 !---------------------------------------------- 422 377 snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass … … 470 425 471 426 !------------------------------------ 472 ! 6) store fields at before time-step427 ! 4) store fields at before time-step 473 428 !------------------------------------ 474 429 ! it is only necessary for the 1st interpolation by Agrif … … 508 463 ! 509 464 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_smi465 TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd 511 466 TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read 512 467 ! 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 468 NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, & 469 & rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, & 470 & rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, & 471 & rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, & 472 & rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, & 473 & sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir 517 474 !!----------------------------------------------------------------------------- 518 475 ! 519 476 REWIND( numnam_ice_ref ) ! Namelist namini in reference namelist : Ice initial state 520 477 READ ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901) 521 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' , lwp)478 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' ) 522 479 REWIND( numnam_ice_cfg ) ! Namelist namini in configuration namelist : Ice initial state 523 480 READ ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 ) 524 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' , lwp)481 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' ) 525 482 IF(lwm) WRITE ( numoni, namini ) 526 483 ! 527 484 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 485 slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi 486 slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms 487 slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd 530 488 ! 531 489 IF(lwp) THEN ! control print … … 534 492 WRITE(numout,*) '~~~~~~~~~~~~~~~' 535 493 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 494 WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini 495 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file 496 WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst 497 IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN 498 WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s 499 WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s 500 WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s 501 WRITE(numout,*) ' initial ice salinity in the north-south rn_smi_ini = ', rn_smi_ini_n,rn_smi_ini_s 502 WRITE(numout,*) ' initial surf temperat in the north-south rn_tsu_ini = ', rn_tsu_ini_n,rn_tsu_ini_s 503 WRITE(numout,*) ' initial ice temperat in the north-south rn_tmi_ini = ', rn_tmi_ini_n,rn_tmi_ini_s 504 WRITE(numout,*) ' initial snw temperat in the north-south rn_tms_ini = ', rn_tms_ini_n,rn_tms_ini_s 505 WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s 506 WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s 507 ENDIF 549 508 ENDIF 550 509 ! … … 554 513 ALLOCATE( si(jpfldi), STAT=ierror ) 555 514 IF( ierror > 0 ) THEN 556 CALL ctl_stop( ' Ice_ini in iceistate: unable to allocate si structure' ) ; RETURN515 CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' ) ; RETURN 557 516 ENDIF 558 517 ! 559 518 DO ifpr = 1, jpfldi 560 519 ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) 561 ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )520 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 562 521 END DO 563 522 ! 564 523 ! 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 ! 524 CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' ) 525 ! 526 ENDIF 527 ! 528 IF( .NOT.ln_pnd ) THEN 529 rn_apd_ini_n = 0. ; rn_apd_ini_s = 0. 530 rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0. 531 CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' ) 569 532 ENDIF 570 533 !
Note: See TracChangeset
for help on using the changeset viewer.
