[821] | 1 | MODULE limthd_2 |
---|
[3] | 2 | !!====================================================================== |
---|
[821] | 3 | !! *** MODULE limthd_2 *** |
---|
[3] | 4 | !! LIM thermo ice model : ice thermodynamic |
---|
| 5 | !!====================================================================== |
---|
[888] | 6 | !! History : 1.0 ! 00-01 (LIM) |
---|
| 7 | !! 2.0 ! 02-07 (C. Ethe, G. Madec) F90 |
---|
| 8 | !! 2.0 ! 03-08 (C. Ethe) add lim_thd_init |
---|
[1218] | 9 | !! - ! 08-2008 (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface |
---|
[888] | 10 | !!--------------------------------------------------------------------- |
---|
[821] | 11 | #if defined key_lim2 |
---|
[3] | 12 | !!---------------------------------------------------------------------- |
---|
[821] | 13 | !! 'key_lim2' : LIM 2.0 sea-ice model |
---|
[3] | 14 | !!---------------------------------------------------------------------- |
---|
[821] | 15 | !! lim_thd_2 : thermodynamic of sea ice |
---|
| 16 | !! lim_thd_init_2 : initialisation of sea-ice thermodynamic |
---|
[3] | 17 | !!---------------------------------------------------------------------- |
---|
[717] | 18 | USE phycst ! physical constants |
---|
[3] | 19 | USE dom_oce ! ocean space and time domain variables |
---|
[717] | 20 | USE lbclnk |
---|
| 21 | USE in_out_manager ! I/O manager |
---|
[821] | 22 | USE ice_2 ! LIM sea-ice variables |
---|
[3] | 23 | USE ice_oce ! sea-ice/ocean variables |
---|
[888] | 24 | USE sbc_oce ! |
---|
| 25 | USE sbc_ice ! |
---|
[821] | 26 | USE thd_ice_2 ! LIM thermodynamic sea-ice variables |
---|
| 27 | USE dom_ice_2 ! LIM sea-ice domain |
---|
| 28 | USE limthd_zdf_2 |
---|
| 29 | USE limthd_lac_2 |
---|
| 30 | USE limtab_2 |
---|
[258] | 31 | USE prtctl ! Print control |
---|
[1218] | 32 | USE cpl_oasis3, ONLY : lk_cpl |
---|
[3] | 33 | |
---|
| 34 | IMPLICIT NONE |
---|
| 35 | PRIVATE |
---|
| 36 | |
---|
[888] | 37 | PUBLIC lim_thd_2 ! called by lim_step |
---|
[3] | 38 | |
---|
[1218] | 39 | REAL(wp) :: epsi20 = 1.e-20 ! constant values |
---|
| 40 | REAL(wp) :: epsi16 = 1.e-16 ! |
---|
| 41 | REAL(wp) :: epsi04 = 1.e-04 ! |
---|
| 42 | REAL(wp) :: rzero = 0.e0 ! |
---|
| 43 | REAL(wp) :: rone = 1.e0 ! |
---|
[70] | 44 | |
---|
[3] | 45 | !! * Substitutions |
---|
[70] | 46 | # include "domzgr_substitute.h90" |
---|
[3] | 47 | # include "vectopt_loop_substitute.h90" |
---|
| 48 | !!-------- ------------------------------------------------------------- |
---|
[1218] | 49 | !! NEMO/LIM 2.0, UCL-LOCEAN-IPSL (2008) |
---|
[1156] | 50 | !! $Id$ |
---|
[888] | 51 | !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) |
---|
[3] | 52 | !!---------------------------------------------------------------------- |
---|
| 53 | |
---|
| 54 | CONTAINS |
---|
| 55 | |
---|
[821] | 56 | SUBROUTINE lim_thd_2( kt ) |
---|
[3] | 57 | !!------------------------------------------------------------------- |
---|
[821] | 58 | !! *** ROUTINE lim_thd_2 *** |
---|
[3] | 59 | !! |
---|
| 60 | !! ** Purpose : This routine manages the ice thermodynamic. |
---|
| 61 | !! |
---|
| 62 | !! ** Action : - Initialisation of some variables |
---|
| 63 | !! - Some preliminary computation (oceanic heat flux |
---|
| 64 | !! at the ice base, snow acc.,heat budget of the leads) |
---|
| 65 | !! - selection of the icy points and put them in an array |
---|
| 66 | !! - call lim_vert_ther for vert ice thermodynamic |
---|
| 67 | !! - back to the geographic grid |
---|
| 68 | !! - selection of points for lateral accretion |
---|
| 69 | !! - call lim_lat_acc for the ice accretion |
---|
| 70 | !! - back to the geographic grid |
---|
| 71 | !! |
---|
[888] | 72 | !! References : Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90 |
---|
[3] | 73 | !!--------------------------------------------------------------------- |
---|
[508] | 74 | INTEGER, INTENT(in) :: kt ! number of iteration |
---|
[888] | 75 | !! |
---|
[1218] | 76 | INTEGER :: ji, jj ! dummy loop indices |
---|
| 77 | INTEGER :: nbpb ! nb of icy pts for thermo. cal. |
---|
| 78 | INTEGER :: nbpac ! nb of pts for lateral accretion |
---|
[258] | 79 | CHARACTER (len=22) :: charout |
---|
[1218] | 80 | REAL(wp) :: zfric_umin = 5e-03 ! lower bound for the friction velocity |
---|
| 81 | REAL(wp) :: zfric_umax = 2e-02 ! upper bound for the friction velocity |
---|
| 82 | REAL(wp) :: zinda ! switch for test. the val. of concen. |
---|
| 83 | REAL(wp) :: zindb, zindg ! switches for test. the val of arg |
---|
| 84 | REAL(wp) :: za , zh, zthsnice ! |
---|
| 85 | REAL(wp) :: zfric_u ! friction velocity |
---|
| 86 | REAL(wp) :: zfnsol ! total non solar heat |
---|
| 87 | REAL(wp) :: zfontn ! heat flux from snow thickness |
---|
| 88 | REAL(wp) :: zfntlat, zpareff ! test. the val. of lead heat budget |
---|
| 89 | REAL(wp) :: zfi ! temporary scalar |
---|
| 90 | REAL(wp), DIMENSION(jpi,jpj) :: zhicifp ! ice thickness for outputs |
---|
| 91 | REAL(wp), DIMENSION(jpi,jpj) :: zqlbsbq ! link with lead energy budget qldif |
---|
[888] | 92 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmsk ! working array |
---|
[3] | 93 | !!------------------------------------------------------------------- |
---|
| 94 | |
---|
[888] | 95 | IF( kt == nit000 ) CALL lim_thd_init_2 ! Initialization (first time-step only) |
---|
[3] | 96 | |
---|
| 97 | !-------------------------------------------! |
---|
| 98 | ! Initilization of diagnostic variables ! |
---|
| 99 | !-------------------------------------------! |
---|
| 100 | |
---|
[1218] | 101 | !!gm needed? yes at least for some of these arrays |
---|
[70] | 102 | rdvosif(:,:) = 0.e0 ! variation of ice volume at surface |
---|
| 103 | rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom |
---|
| 104 | fdvolif(:,:) = 0.e0 ! total variation of ice volume |
---|
| 105 | rdvonif(:,:) = 0.e0 ! lateral variation of ice volume |
---|
| 106 | fstric (:,:) = 0.e0 ! part of solar radiation absorbing inside the ice |
---|
| 107 | fscmbq (:,:) = 0.e0 ! linked with fstric |
---|
| 108 | ffltbif(:,:) = 0.e0 ! linked with fstric |
---|
| 109 | qfvbq (:,:) = 0.e0 ! linked with fstric |
---|
| 110 | rdmsnif(:,:) = 0.e0 ! variation of snow mass per unit area |
---|
| 111 | rdmicif(:,:) = 0.e0 ! variation of ice mass per unit area |
---|
| 112 | hicifp (:,:) = 0.e0 ! daily thermodynamic ice production. |
---|
[258] | 113 | zmsk (:,:,:) = 0.e0 |
---|
[3] | 114 | |
---|
[1218] | 115 | ! set to zero snow thickness smaller than epsi04 |
---|
[3] | 116 | DO jj = 1, jpj |
---|
| 117 | DO ji = 1, jpi |
---|
[1218] | 118 | hsnif(ji,jj) = hsnif(ji,jj) * MAX( rzero, SIGN( rone , hsnif(ji,jj) - epsi04 ) ) |
---|
[3] | 119 | END DO |
---|
| 120 | END DO |
---|
[1218] | 121 | !!gm better coded (do not use SIGN...) |
---|
| 122 | ! WHERE( hsnif(:,:) < epsi04 ) hsnif(:,:) = 0.e0 |
---|
| 123 | !!gm |
---|
[258] | 124 | |
---|
[1218] | 125 | IF(ln_ctl) CALL prt_ctl( tab2d_1=hsnif, clinfo1=' lim_thd: hsnif : ' ) |
---|
[3] | 126 | |
---|
| 127 | !-----------------------------------! |
---|
| 128 | ! Treatment of particular cases ! |
---|
| 129 | !-----------------------------------! |
---|
| 130 | |
---|
| 131 | DO jj = 1, jpj |
---|
| 132 | DO ji = 1, jpi |
---|
| 133 | ! snow is transformed into ice if the original ice cover disappears. |
---|
[1218] | 134 | zindg = tms(ji,jj) * MAX( rzero , SIGN( rone , -hicif(ji,jj) ) ) |
---|
[3] | 135 | hicif(ji,jj) = hicif(ji,jj) + zindg * rhosn * hsnif(ji,jj) / rau0 |
---|
[1218] | 136 | hsnif(ji,jj) = ( rone - zindg ) * hsnif(ji,jj) + zindg * hicif(ji,jj) * ( rau0 - rhoic ) / rhosn |
---|
[3] | 137 | dmgwi(ji,jj) = zindg * (1.0 - frld(ji,jj)) * rhoic * hicif(ji,jj) ! snow/ice mass |
---|
| 138 | |
---|
| 139 | ! the lead fraction, frld, must be little than or equal to amax (ice ridging). |
---|
| 140 | zthsnice = hsnif(ji,jj) + hicif(ji,jj) |
---|
[1218] | 141 | zindb = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) ) |
---|
| 142 | za = zindb * MIN( rone, ( 1.0 - frld(ji,jj) ) * uscomi ) |
---|
[3] | 143 | hsnif (ji,jj) = hsnif(ji,jj) * za |
---|
| 144 | hicif (ji,jj) = hicif(ji,jj) * za |
---|
| 145 | qstoif(ji,jj) = qstoif(ji,jj) * za |
---|
[1218] | 146 | frld (ji,jj) = 1.0 - zindb * ( 1.0 - frld(ji,jj) ) / MAX( za, epsi20 ) |
---|
[3] | 147 | |
---|
| 148 | ! the in situ ice thickness, hicif, must be equal to or greater than hiclim. |
---|
[1218] | 149 | zh = MAX( rone , zindb * hiclim / MAX( hicif(ji,jj), epsi20 ) ) |
---|
[3] | 150 | hsnif (ji,jj) = hsnif(ji,jj) * zh |
---|
| 151 | hicif (ji,jj) = hicif(ji,jj) * zh |
---|
| 152 | qstoif(ji,jj) = qstoif(ji,jj) * zh |
---|
| 153 | frld (ji,jj) = ( frld(ji,jj) + ( zh - 1.0 ) ) / zh |
---|
| 154 | END DO |
---|
| 155 | END DO |
---|
[258] | 156 | |
---|
| 157 | IF(ln_ctl) THEN |
---|
[1218] | 158 | CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif : ' ) |
---|
| 159 | CALL prt_ctl( tab2d_1=hsnif , clinfo1=' lim_thd: hsnif : ' ) |
---|
| 160 | CALL prt_ctl( tab2d_1=dmgwi , clinfo1=' lim_thd: dmgwi : ' ) |
---|
| 161 | CALL prt_ctl( tab2d_1=qstoif, clinfo1=' lim_thd: qstoif : ' ) |
---|
| 162 | CALL prt_ctl( tab2d_1=frld , clinfo1=' lim_thd: frld : ' ) |
---|
[3] | 163 | ENDIF |
---|
| 164 | |
---|
| 165 | |
---|
| 166 | !-------------------------------! |
---|
| 167 | ! Thermodynamics of sea ice ! |
---|
| 168 | !-------------------------------! |
---|
| 169 | |
---|
| 170 | ! Partial computation of forcing for the thermodynamic sea ice model. |
---|
| 171 | !-------------------------------------------------------------------------- |
---|
| 172 | |
---|
[888] | 173 | sst_m(:,:) = sst_m(:,:) + rt0 |
---|
| 174 | |
---|
[3] | 175 | !CDIR NOVERRCHK |
---|
| 176 | DO jj = 1, jpj |
---|
| 177 | !CDIR NOVERRCHK |
---|
| 178 | DO ji = 1, jpi |
---|
| 179 | zthsnice = hsnif(ji,jj) + hicif(ji,jj) |
---|
[1218] | 180 | zindb = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) ) |
---|
[3] | 181 | pfrld(ji,jj) = frld(ji,jj) |
---|
[1218] | 182 | zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) |
---|
[3] | 183 | |
---|
| 184 | ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget |
---|
[70] | 185 | thcm(ji,jj) = 0.e0 |
---|
[3] | 186 | |
---|
| 187 | ! net downward heat flux from the ice to the ocean, expressed as a function of ocean |
---|
| 188 | ! temperature and turbulent mixing (McPhee, 1992) |
---|
| 189 | zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) ! friction velocity |
---|
[888] | 190 | fdtcn(ji,jj) = zindb * rau0 * rcp * 0.006 * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) ) |
---|
[3] | 191 | qdtcn(ji,jj) = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice |
---|
| 192 | |
---|
| 193 | ! partial computation of the lead energy budget (qldif) |
---|
[1218] | 194 | #if defined key_coupled |
---|
| 195 | zfi = 1.0 - pfrld(ji,jj) |
---|
| 196 | qldif(ji,jj) = tms(ji,jj) * rdt_ice & |
---|
[1463] | 197 | & * ( ( qsr_tot(ji,jj) - qsr_ice(ji,jj,1) * zfi ) * ( 1.0 - thcm(ji,jj) ) & |
---|
| 198 | & + ( qns_tot(ji,jj) - qns_ice(ji,jj,1) * zfi ) & |
---|
[1218] | 199 | & + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) ) ) |
---|
| 200 | #else |
---|
| 201 | zfontn = ( sprecip(ji,jj) / rhosn ) * xlsn ! energy for melting solid precipitation |
---|
[888] | 202 | zfnsol = qns(ji,jj) ! total non solar flux over the ocean |
---|
| 203 | qldif(ji,jj) = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) ) & |
---|
[3] | 204 | & + zfnsol + fdtcn(ji,jj) - zfontn & |
---|
| 205 | & + ( 1.0 - zindb ) * fsbbq(ji,jj) ) & |
---|
[1218] | 206 | & * frld(ji,jj) * rdt_ice |
---|
| 207 | !!$ qldif(ji,jj) = tms(ji,jj) * rdt_ice * frld(ji,jj) |
---|
| 208 | !!$ & * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) ) & |
---|
| 209 | !!$ & + qns(ji,jj) + fdtcn(ji,jj) - zfontn & |
---|
| 210 | !!$ & + ( 1.0 - zindb ) * fsbbq(ji,jj) ) & |
---|
| 211 | #endif |
---|
[3] | 212 | ! parlat : percentage of energy used for lateral ablation (0.0) |
---|
[1218] | 213 | zfntlat = 1.0 - MAX( rzero , SIGN( rone , - qldif(ji,jj) ) ) |
---|
[3] | 214 | zpareff = 1.0 + ( parlat - 1.0 ) * zinda * zfntlat |
---|
| 215 | zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( (1.0 - frld(ji,jj)) * rdt_ice , epsi16 ) |
---|
| 216 | qldif (ji,jj) = zpareff * qldif(ji,jj) |
---|
| 217 | qdtcn (ji,jj) = zpareff * qdtcn(ji,jj) |
---|
| 218 | |
---|
| 219 | ! energy needed to bring ocean surface layer until its freezing |
---|
[888] | 220 | qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda ) |
---|
[3] | 221 | |
---|
| 222 | ! calculate oceanic heat flux. |
---|
| 223 | fbif (ji,jj) = zindb * ( fsbbq(ji,jj) / MAX( (1.0 - frld(ji,jj)) , epsi20 ) + fdtcn(ji,jj) ) |
---|
| 224 | |
---|
| 225 | ! computation of the daily thermodynamic ice production (only needed for output) |
---|
| 226 | zhicifp(ji,jj) = hicif(ji,jj) * ( 1.0 - frld(ji,jj) ) |
---|
| 227 | END DO |
---|
| 228 | END DO |
---|
| 229 | |
---|
[888] | 230 | sst_m(:,:) = sst_m(:,:) - rt0 |
---|
[3] | 231 | |
---|
| 232 | ! Select icy points and fulfill arrays for the vectorial grid. |
---|
| 233 | !---------------------------------------------------------------------- |
---|
| 234 | nbpb = 0 |
---|
| 235 | DO jj = 1, jpj |
---|
| 236 | DO ji = 1, jpi |
---|
| 237 | IF ( frld(ji,jj) < 1.0 ) THEN |
---|
| 238 | nbpb = nbpb + 1 |
---|
| 239 | npb(nbpb) = (jj - 1) * jpi + ji |
---|
| 240 | ENDIF |
---|
| 241 | END DO |
---|
| 242 | END DO |
---|
[258] | 243 | |
---|
| 244 | IF(ln_ctl) THEN |
---|
| 245 | CALL prt_ctl(tab2d_1=pfrld, clinfo1=' lim_thd: pfrld : ', tab2d_2=thcm , clinfo2=' thcm : ') |
---|
| 246 | CALL prt_ctl(tab2d_1=fdtcn, clinfo1=' lim_thd: fdtcn : ', tab2d_2=qdtcn , clinfo2=' qdtcn : ') |
---|
| 247 | CALL prt_ctl(tab2d_1=qldif, clinfo1=' lim_thd: qldif : ', tab2d_2=zqlbsbq, clinfo2=' zqlbsbq : ') |
---|
| 248 | CALL prt_ctl(tab2d_1=qcmif, clinfo1=' lim_thd: qcmif : ', tab2d_2=fbif , clinfo2=' fbif : ') |
---|
| 249 | zmsk(:,:,1) = tms(:,:) |
---|
| 250 | CALL prt_ctl(tab2d_1=qcmif , clinfo1=' lim_thd: qcmif : ', mask1=zmsk) |
---|
| 251 | CALL prt_ctl(tab2d_1=zhicifp, clinfo1=' lim_thd: zhicifp : ') |
---|
| 252 | WRITE(charout, FMT="('lim_thd: nbpb = ',I4)") nbpb |
---|
| 253 | CALL prt_ctl_info(charout) |
---|
[3] | 254 | ENDIF |
---|
| 255 | |
---|
| 256 | |
---|
| 257 | ! If there is no ice, do nothing. Otherwise, compute Top and Bottom accretion/ablation |
---|
| 258 | !------------------------------------------------------------------------------------ |
---|
| 259 | |
---|
[1218] | 260 | IF( nbpb > 0 ) THEN |
---|
| 261 | ! |
---|
[3] | 262 | ! put the variable in a 1-D array for thermodynamics process |
---|
[1463] | 263 | CALL tab_2d_1d_2( nbpb, frld_1d (1:nbpb) , frld , jpi, jpj, npb(1:nbpb) ) |
---|
| 264 | CALL tab_2d_1d_2( nbpb, h_ice_1d (1:nbpb) , hicif , jpi, jpj, npb(1:nbpb) ) |
---|
| 265 | CALL tab_2d_1d_2( nbpb, h_snow_1d (1:nbpb) , hsnif , jpi, jpj, npb(1:nbpb) ) |
---|
| 266 | CALL tab_2d_1d_2( nbpb, sist_1d (1:nbpb) , sist , jpi, jpj, npb(1:nbpb) ) |
---|
| 267 | CALL tab_2d_1d_2( nbpb, tbif_1d (1:nbpb , 1 ), tbif(:,:,1) , jpi, jpj, npb(1:nbpb) ) |
---|
| 268 | CALL tab_2d_1d_2( nbpb, tbif_1d (1:nbpb , 2 ), tbif(:,:,2) , jpi, jpj, npb(1:nbpb) ) |
---|
| 269 | CALL tab_2d_1d_2( nbpb, tbif_1d (1:nbpb , 3 ), tbif(:,:,3) , jpi, jpj, npb(1:nbpb) ) |
---|
| 270 | CALL tab_2d_1d_2( nbpb, qsr_ice_1d (1:nbpb) , qsr_ice(:,:,1) , jpi, jpj, npb(1:nbpb) ) |
---|
| 271 | CALL tab_2d_1d_2( nbpb, fr1_i0_1d (1:nbpb) , fr1_i0 , jpi, jpj, npb(1:nbpb) ) |
---|
| 272 | CALL tab_2d_1d_2( nbpb, fr2_i0_1d (1:nbpb) , fr2_i0 , jpi, jpj, npb(1:nbpb) ) |
---|
| 273 | CALL tab_2d_1d_2( nbpb, qns_ice_1d(1:nbpb) , qns_ice(:,:,1), jpi, jpj, npb(1:nbpb) ) |
---|
| 274 | CALL tab_2d_1d_2( nbpb, dqns_ice_1d(1:nbpb) , dqns_ice(:,:,1), jpi, jpj, npb(1:nbpb) ) |
---|
[1218] | 275 | IF( .NOT. lk_cpl ) THEN |
---|
[1463] | 276 | CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb) , qla_ice(:,:,1), jpi, jpj, npb(1:nbpb) ) |
---|
| 277 | CALL tab_2d_1d_2( nbpb, dqla_ice_1d(1:nbpb) , dqla_ice(:,:,1), jpi, jpj, npb(1:nbpb) ) |
---|
[1218] | 278 | ENDIF |
---|
[821] | 279 | CALL tab_2d_1d_2( nbpb, tfu_1d (1:nbpb) , tfu , jpi, jpj, npb(1:nbpb) ) |
---|
| 280 | CALL tab_2d_1d_2( nbpb, sprecip_1d (1:nbpb) , sprecip , jpi, jpj, npb(1:nbpb) ) |
---|
| 281 | CALL tab_2d_1d_2( nbpb, fbif_1d (1:nbpb) , fbif , jpi, jpj, npb(1:nbpb) ) |
---|
| 282 | CALL tab_2d_1d_2( nbpb, thcm_1d (1:nbpb) , thcm , jpi, jpj, npb(1:nbpb) ) |
---|
| 283 | CALL tab_2d_1d_2( nbpb, qldif_1d (1:nbpb) , qldif , jpi, jpj, npb(1:nbpb) ) |
---|
| 284 | CALL tab_2d_1d_2( nbpb, qstbif_1d (1:nbpb) , qstoif , jpi, jpj, npb(1:nbpb) ) |
---|
| 285 | CALL tab_2d_1d_2( nbpb, rdmicif_1d (1:nbpb) , rdmicif , jpi, jpj, npb(1:nbpb) ) |
---|
| 286 | CALL tab_2d_1d_2( nbpb, dmgwi_1d (1:nbpb) , dmgwi , jpi, jpj, npb(1:nbpb) ) |
---|
| 287 | CALL tab_2d_1d_2( nbpb, qlbbq_1d (1:nbpb) , zqlbsbq , jpi, jpj, npb(1:nbpb) ) |
---|
[1218] | 288 | ! |
---|
[821] | 289 | CALL lim_thd_zdf_2( 1, nbpb ) ! compute ice growth |
---|
[1218] | 290 | ! |
---|
[3] | 291 | ! back to the geographic grid. |
---|
[821] | 292 | CALL tab_1d_2d_2( nbpb, frld , npb, frld_1d (1:nbpb) , jpi, jpj ) |
---|
| 293 | CALL tab_1d_2d_2( nbpb, hicif , npb, h_ice_1d (1:nbpb) , jpi, jpj ) |
---|
| 294 | CALL tab_1d_2d_2( nbpb, hsnif , npb, h_snow_1d (1:nbpb) , jpi, jpj ) |
---|
| 295 | CALL tab_1d_2d_2( nbpb, sist , npb, sist_1d (1:nbpb) , jpi, jpj ) |
---|
| 296 | CALL tab_1d_2d_2( nbpb, tbif(:,:,1), npb, tbif_1d (1:nbpb , 1 ), jpi, jpj ) |
---|
| 297 | CALL tab_1d_2d_2( nbpb, tbif(:,:,2), npb, tbif_1d (1:nbpb , 2 ), jpi, jpj ) |
---|
| 298 | CALL tab_1d_2d_2( nbpb, tbif(:,:,3), npb, tbif_1d (1:nbpb , 3 ), jpi, jpj ) |
---|
| 299 | CALL tab_1d_2d_2( nbpb, fscmbq , npb, fscbq_1d (1:nbpb) , jpi, jpj ) |
---|
| 300 | CALL tab_1d_2d_2( nbpb, ffltbif , npb, fltbif_1d (1:nbpb) , jpi, jpj ) |
---|
| 301 | CALL tab_1d_2d_2( nbpb, fstric , npb, fstbif_1d (1:nbpb) , jpi, jpj ) |
---|
| 302 | CALL tab_1d_2d_2( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj ) |
---|
| 303 | CALL tab_1d_2d_2( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj ) |
---|
| 304 | CALL tab_1d_2d_2( nbpb, qstoif , npb, qstbif_1d (1:nbpb) , jpi, jpj ) |
---|
| 305 | CALL tab_1d_2d_2( nbpb, rdmicif , npb, rdmicif_1d(1:nbpb) , jpi, jpj ) |
---|
| 306 | CALL tab_1d_2d_2( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj ) |
---|
| 307 | CALL tab_1d_2d_2( nbpb, rdmsnif , npb, rdmsnif_1d(1:nbpb) , jpi, jpj ) |
---|
| 308 | CALL tab_1d_2d_2( nbpb, rdvosif , npb, dvsbq_1d (1:nbpb) , jpi, jpj ) |
---|
| 309 | CALL tab_1d_2d_2( nbpb, rdvobif , npb, dvbbq_1d (1:nbpb) , jpi, jpj ) |
---|
| 310 | CALL tab_1d_2d_2( nbpb, fdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj ) |
---|
| 311 | CALL tab_1d_2d_2( nbpb, rdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj ) |
---|
[1218] | 312 | ! |
---|
[3] | 313 | ENDIF |
---|
| 314 | |
---|
| 315 | |
---|
[1218] | 316 | ! Up-date sea ice thickness |
---|
| 317 | !-------------------------- |
---|
[3] | 318 | DO jj = 1, jpj |
---|
| 319 | DO ji = 1, jpi |
---|
| 320 | phicif(ji,jj) = hicif(ji,jj) |
---|
[1218] | 321 | hicif(ji,jj) = hicif(ji,jj) * ( rone - MAX( rzero, SIGN( rone, - ( 1.0 - frld(ji,jj) ) ) ) ) |
---|
[3] | 322 | END DO |
---|
| 323 | END DO |
---|
| 324 | |
---|
| 325 | |
---|
[1218] | 326 | ! Tricky trick : add 2 to frld in the Southern Hemisphere |
---|
| 327 | !-------------------------------------------------------- |
---|
[192] | 328 | IF( fcor(1,1) < 0.e0 ) THEN |
---|
[421] | 329 | DO jj = 1, njeqm1 |
---|
[192] | 330 | DO ji = 1, jpi |
---|
| 331 | frld(ji,jj) = frld(ji,jj) + 2.0 |
---|
| 332 | END DO |
---|
[3] | 333 | END DO |
---|
[192] | 334 | ENDIF |
---|
[3] | 335 | |
---|
| 336 | |
---|
[1218] | 337 | ! Select points for lateral accretion (this occurs when heat exchange |
---|
| 338 | ! between ice and ocean is negative; ocean losing heat) |
---|
[3] | 339 | !----------------------------------------------------------------- |
---|
| 340 | nbpac = 0 |
---|
| 341 | DO jj = 1, jpj |
---|
| 342 | DO ji = 1, jpi |
---|
[70] | 343 | !i yes! IF ( ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN |
---|
| 344 | IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN |
---|
[3] | 345 | nbpac = nbpac + 1 |
---|
| 346 | npac( nbpac ) = (jj - 1) * jpi + ji |
---|
| 347 | ENDIF |
---|
| 348 | END DO |
---|
| 349 | END DO |
---|
| 350 | |
---|
[258] | 351 | IF(ln_ctl) THEN |
---|
| 352 | CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif : ', tab2d_2=hicif, clinfo2=' hicif : ') |
---|
| 353 | WRITE(charout, FMT="('lim_thd: nbpac = ',I4)") nbpac |
---|
| 354 | CALL prt_ctl_info(charout) |
---|
[3] | 355 | ENDIF |
---|
| 356 | |
---|
[1218] | 357 | |
---|
| 358 | ! If ocean gains heat do nothing ; otherwise, one performs lateral accretion |
---|
[3] | 359 | !-------------------------------------------------------------------------------- |
---|
[70] | 360 | IF( nbpac > 0 ) THEN |
---|
[1218] | 361 | ! |
---|
[3] | 362 | !...Put the variable in a 1-D array for lateral accretion |
---|
[821] | 363 | CALL tab_2d_1d_2( nbpac, frld_1d (1:nbpac) , frld , jpi, jpj, npac(1:nbpac) ) |
---|
| 364 | CALL tab_2d_1d_2( nbpac, h_snow_1d (1:nbpac) , hsnif , jpi, jpj, npac(1:nbpac) ) |
---|
| 365 | CALL tab_2d_1d_2( nbpac, h_ice_1d (1:nbpac) , hicif , jpi, jpj, npac(1:nbpac) ) |
---|
| 366 | CALL tab_2d_1d_2( nbpac, tbif_1d (1:nbpac , 1 ), tbif(:,:,1), jpi, jpj, npac(1:nbpac) ) |
---|
| 367 | CALL tab_2d_1d_2( nbpac, tbif_1d (1:nbpac , 2 ), tbif(:,:,2), jpi, jpj, npac(1:nbpac) ) |
---|
| 368 | CALL tab_2d_1d_2( nbpac, tbif_1d (1:nbpac , 3 ), tbif(:,:,3), jpi, jpj, npac(1:nbpac) ) |
---|
| 369 | CALL tab_2d_1d_2( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) ) |
---|
| 370 | CALL tab_2d_1d_2( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) ) |
---|
| 371 | CALL tab_2d_1d_2( nbpac, qstbif_1d (1:nbpac) , qstoif , jpi, jpj, npac(1:nbpac) ) |
---|
| 372 | CALL tab_2d_1d_2( nbpac, rdmicif_1d(1:nbpac) , rdmicif , jpi, jpj, npac(1:nbpac) ) |
---|
| 373 | CALL tab_2d_1d_2( nbpac, dvlbq_1d (1:nbpac) , fdvolif , jpi, jpj, npac(1:nbpac) ) |
---|
| 374 | CALL tab_2d_1d_2( nbpac, tfu_1d (1:nbpac) , tfu , jpi, jpj, npac(1:nbpac) ) |
---|
[1218] | 375 | ! |
---|
| 376 | CALL lim_thd_lac_2( 1 , nbpac ) ! lateral accretion routine. |
---|
| 377 | ! |
---|
[70] | 378 | ! back to the geographic grid |
---|
[821] | 379 | CALL tab_1d_2d_2( nbpac, frld , npac(1:nbpac), frld_1d (1:nbpac) , jpi, jpj ) |
---|
| 380 | CALL tab_1d_2d_2( nbpac, hsnif , npac(1:nbpac), h_snow_1d (1:nbpac) , jpi, jpj ) |
---|
| 381 | CALL tab_1d_2d_2( nbpac, hicif , npac(1:nbpac), h_ice_1d (1:nbpac) , jpi, jpj ) |
---|
| 382 | CALL tab_1d_2d_2( nbpac, tbif(:,:,1), npac(1:nbpac), tbif_1d (1:nbpac , 1 ), jpi, jpj ) |
---|
| 383 | CALL tab_1d_2d_2( nbpac, tbif(:,:,2), npac(1:nbpac), tbif_1d (1:nbpac , 2 ), jpi, jpj ) |
---|
| 384 | CALL tab_1d_2d_2( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d (1:nbpac , 3 ), jpi, jpj ) |
---|
| 385 | CALL tab_1d_2d_2( nbpac, qstoif , npac(1:nbpac), qstbif_1d (1:nbpac) , jpi, jpj ) |
---|
| 386 | CALL tab_1d_2d_2( nbpac, rdmicif , npac(1:nbpac), rdmicif_1d(1:nbpac) , jpi, jpj ) |
---|
| 387 | CALL tab_1d_2d_2( nbpac, fdvolif , npac(1:nbpac), dvlbq_1d (1:nbpac) , jpi, jpj ) |
---|
[1218] | 388 | ! |
---|
[70] | 389 | ENDIF |
---|
[3] | 390 | |
---|
| 391 | |
---|
[1218] | 392 | ! Recover frld values between 0 and 1 in the Southern Hemisphere (tricky trick) |
---|
| 393 | ! Update daily thermodynamic ice production. |
---|
[70] | 394 | !------------------------------------------------------------------------------ |
---|
[3] | 395 | DO jj = 1, jpj |
---|
| 396 | DO ji = 1, jpi |
---|
| 397 | frld (ji,jj) = MIN( frld(ji,jj), ABS( frld(ji,jj) - 2.0 ) ) |
---|
| 398 | hicifp(ji,jj) = hicif(ji,jj) * ( 1.0 - frld(ji,jj) ) - zhicifp(ji,jj) + hicifp(ji,jj) |
---|
| 399 | END DO |
---|
| 400 | END DO |
---|
| 401 | |
---|
[258] | 402 | IF(ln_ctl) THEN |
---|
| 403 | CALL prt_ctl_info(' lim_thd end ') |
---|
[1218] | 404 | CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif : ', tab2d_2=hsnif , clinfo2=' hsnif : ' ) |
---|
| 405 | CALL prt_ctl( tab2d_1=frld , clinfo1=' lim_thd: frld : ', tab2d_2=hicifp, clinfo2=' hicifp : ' ) |
---|
| 406 | CALL prt_ctl( tab2d_1=phicif , clinfo1=' lim_thd: phicif : ', tab2d_2=pfrld , clinfo2=' pfrld : ' ) |
---|
| 407 | CALL prt_ctl( tab2d_1=sist , clinfo1=' lim_thd: sist : ' ) |
---|
| 408 | CALL prt_ctl( tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1 : ' ) |
---|
| 409 | CALL prt_ctl( tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2 : ' ) |
---|
| 410 | CALL prt_ctl( tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3 : ' ) |
---|
| 411 | CALL prt_ctl( tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn : ', tab2d_2=qdtcn , clinfo2=' qdtcn : ' ) |
---|
| 412 | CALL prt_ctl( tab2d_1=qstoif , clinfo1=' lim_thd: qstoif : ', tab2d_2=fsbbq , clinfo2=' fsbbq : ' ) |
---|
[3] | 413 | ENDIF |
---|
[888] | 414 | ! |
---|
[821] | 415 | END SUBROUTINE lim_thd_2 |
---|
[3] | 416 | |
---|
[719] | 417 | |
---|
[821] | 418 | SUBROUTINE lim_thd_init_2 |
---|
[3] | 419 | !!------------------------------------------------------------------- |
---|
[821] | 420 | !! *** ROUTINE lim_thd_init_2 *** |
---|
[3] | 421 | !! |
---|
| 422 | !! ** Purpose : Physical constants and parameters linked to the ice |
---|
| 423 | !! thermodynamics |
---|
| 424 | !! |
---|
| 425 | !! ** Method : Read the namicethd namelist and check the ice-thermo |
---|
| 426 | !! parameter values called at the first timestep (nit000) |
---|
| 427 | !! |
---|
| 428 | !! ** input : Namelist namicether |
---|
| 429 | !!------------------------------------------------------------------- |
---|
| 430 | NAMELIST/namicethd/ hmelt , hiccrit, hicmin, hiclim, amax , & |
---|
| 431 | & swiqst, sbeta , parlat, hakspl, hibspl, exld, & |
---|
| 432 | & hakdif, hnzst , thth , parsub, alphs |
---|
| 433 | !!------------------------------------------------------------------- |
---|
[1218] | 434 | ! |
---|
| 435 | REWIND( numnam_ice ) ! read namelist |
---|
[3] | 436 | READ ( numnam_ice , namicethd ) |
---|
[1218] | 437 | IF( lk_cpl .AND. parsub /= 0.0 ) CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' ) |
---|
| 438 | ! |
---|
| 439 | IF(lwp) THEN ! control print |
---|
[3] | 440 | WRITE(numout,*) |
---|
[821] | 441 | WRITE(numout,*)'lim_thd_init_2: ice parameters for ice thermodynamic computation ' |
---|
| 442 | WRITE(numout,*)'~~~~~~~~~~~~~~' |
---|
[70] | 443 | WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt |
---|
| 444 | WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2) = ', hiccrit |
---|
| 445 | WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin |
---|
| 446 | WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim |
---|
| 447 | WRITE(numout,*)' maximum lead fraction amax = ', amax |
---|
[1218] | 448 | WRITE(numout,*)' energy stored in brine pocket (=1) or not (=0) swiqst = ', swiqst |
---|
[70] | 449 | WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' |
---|
| 450 | WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta |
---|
| 451 | WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat |
---|
| 452 | WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl |
---|
| 453 | WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl |
---|
| 454 | WRITE(numout,*)' exponent for leads-closure rate exld = ', exld |
---|
| 455 | WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif |
---|
| 456 | WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth |
---|
| 457 | WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst |
---|
| 458 | WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub |
---|
| 459 | WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs |
---|
[3] | 460 | ENDIF |
---|
[1218] | 461 | ! |
---|
[3] | 462 | uscomi = 1.0 / ( 1.0 - amax ) ! inverse of minimum lead fraction |
---|
| 463 | rcdsn = hakdif * rcdsn |
---|
| 464 | rcdic = hakdif * rcdic |
---|
[1218] | 465 | ! |
---|
| 466 | IF( hsndif > 100.e0 .OR. hicdif > 100.e0 ) THEN |
---|
[3] | 467 | cnscg = 0.e0 |
---|
| 468 | ELSE |
---|
| 469 | cnscg = rcpsn / rcpic ! ratio rcpsn/rcpic |
---|
| 470 | ENDIF |
---|
[1218] | 471 | ! |
---|
[821] | 472 | END SUBROUTINE lim_thd_init_2 |
---|
[3] | 473 | |
---|
| 474 | #else |
---|
[70] | 475 | !!---------------------------------------------------------------------- |
---|
[821] | 476 | !! Default option Dummy module NO LIM 2.0 sea-ice model |
---|
[70] | 477 | !!---------------------------------------------------------------------- |
---|
[3] | 478 | CONTAINS |
---|
[821] | 479 | SUBROUTINE lim_thd_2 ! Dummy routine |
---|
| 480 | END SUBROUTINE lim_thd_2 |
---|
[3] | 481 | #endif |
---|
| 482 | |
---|
| 483 | !!====================================================================== |
---|
[821] | 484 | END MODULE limthd_2 |
---|