Changeset 4161
- Timestamp:
- 2013-11-07T11:01:27+01:00 (11 years ago)
- Location:
- branches/2013/dev_LOCEAN_2013/NEMOGCM
- Files:
-
- 9 added
- 8 deleted
- 49 edited
- 6 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef_default.xml
r4153 r4161 26 26 27 27 <file_group id="1h" output_freq="1h" output_level="10" enabled=".TRUE."/> <!-- 1h files --> 28 28 29 <file_group id="2h" output_freq="2h" output_level="10" enabled=".TRUE."/> <!-- 2h files --> 30 29 31 <file_group id="3h" output_freq="3h" output_level="10" enabled=".TRUE."/> <!-- 3h files --> 32 30 33 <file_group id="4h" output_freq="4h" output_level="10" enabled=".TRUE."/> <!-- 4h files --> 31 <file_group id="6h" output_freq="6h" output_level="10" enabled=".TRUE."/> <!-- 6h files --> 32 33 <file_group id="1d" output_freq="1d" output_level="10" enabled=".TRUE."> <!-- 1d files --> 34 35 <file_group id="6h" output_freq="6h" output_level="10" enabled=".TRUE."/> <!-- 6h files --> 36 37 <file_group id="1d" output_freq="1d" output_level="10" enabled=".TRUE." > <!-- 1d files --> 34 38 35 39 <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" > … … 37 41 <field field_ref="sss" name="sos" long_name="sea_surface_salinity" /> 38 42 <field field_ref="ssh" name="zos" long_name="sea_surface_height_above_geoid" /> 39 </file>40 41 <file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >42 <field field_ref="suoce" name="uos" long_name="sea_surface_x_velocity" />43 </file>44 45 <file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >46 <field field_ref="svoce" name="vos" long_name="sea_surface_y_velocity" />47 </file>48 49 </file_group>50 51 <file_group id="3d" output_freq="3d" output_level="10" enabled=".TRUE."/> <!-- 3d files -->52 53 <file_group id="5d" output_freq="5d" output_level="10" enabled=".TRUE."> <!-- 5d files -->54 55 <file id="file4" name_suffix="_grid_T" description="ocean T grid variables" >56 43 <field field_ref="toce" name="thetao" long_name="sea_water_potential_temperature" /> 57 44 <field field_ref="soce" name="so" long_name="sea_water_salinity" /> 58 <field field_ref="sst" name="tos" long_name="sea_surface_temperature" />59 45 <field field_ref="sst2" name="tossq" long_name="square_of_sea_surface_temperature" /> 60 <field field_ref="sss" name="sos" long_name="sea_surface_salinity" />61 <field field_ref="ssh" name="zos" long_name="sea_surface_height_above_geoid" />62 46 <field field_ref="ssh2" name="zossq" long_name="square_of_sea_surface_height_above_geoid" /> 47 <field field_ref="mldkz5" /> 48 <field field_ref="mldr10_1" /> 49 </file> 50 51 <file id="file2" name_suffix="_SBC" description="surface fluxes variables" > <!-- time step automaticaly defined based on nn_fsbc --> 63 52 <field field_ref="empmr" name="wfo" long_name="water_flux_into_sea_water" /> 64 53 <field field_ref="qsr" name="rsntds" long_name="surface_net_downward_shortwave_flux" /> 65 54 <field field_ref="qt" name="tohfls" long_name="surface_net_downward_total_heat_flux" /> 66 <field field_ref="taum" /> 67 <field field_ref="mldkz5" /> 68 <field field_ref="mldr10_1" /> 69 </file> 70 71 <file id="file5" name_suffix="_grid_U" description="ocean U grid variables" > 55 <field field_ref="saltflx" name="sosflxdo" /> 56 <field field_ref="taum" name="taum" /> 57 <field field_ref="wspd" name="sowindsp" /> 58 <field field_ref="precip" name="soprecip" /> 59 <!-- ice and snow --> 60 <field field_ref="snowpre" /> 61 <field field_ref="qsr_io" name="iicesflx" /> 62 <field field_ref="qns_io" name="iicenflx" /> 63 <field field_ref="utau_ice" name="iicestru" /> 64 <field field_ref="vtau_ice" name="iicestrv" /> 65 66 </file> 67 68 <file id="file3" name_suffix="_grid_U" description="ocean U grid variables" > 69 <field field_ref="suoce" name="uos" long_name="sea_surface_x_velocity" /> 72 70 <field field_ref="uoce" name="uo" long_name="sea_water_x_velocity" /> 73 <field field_ref="suoce" name="uos" long_name="sea_surface_x_velocity" />74 71 <field field_ref="utau" name="tauuo" long_name="surface_downward_x_stress" /> 75 </file> 76 77 <file id="file6" name_suffix="_grid_V" description="ocean V grid variables" > 72 <!-- available with key_diaar5 --> 73 <field field_ref="u_masstr" name="vozomatr" /> 74 <field field_ref="u_heattr" name="sozohetr" /> 75 </file> 76 77 <file id="file4" name_suffix="_grid_V" description="ocean V grid variables" > 78 <field field_ref="svoce" name="vos" long_name="sea_surface_y_velocity" /> 78 79 <field field_ref="voce" name="vo" long_name="sea_water_y_velocity" /> 79 <field field_ref="svoce" name="vos" long_name="sea_surface_y_velocity" />80 80 <field field_ref="vtau" name="tauvo" long_name="surface_downward_y_stress" /> 81 </file> 82 83 <file id="file7" name_suffix="_grid_W" description="ocean W grid variables" > 81 <!-- available with key_diaar5 --> 82 <field field_ref="v_masstr" name="vomematr" /> 83 <field field_ref="v_heattr" name="somehetr" /> 84 </file> 85 86 <file id="file5" name_suffix="_grid_W" description="ocean W grid variables" > 84 87 <field field_ref="woce" name="wo" long_name="ocean vertical velocity" /> 85 88 <field field_ref="avt" name="difvho" long_name="ocean_vertical_heat_diffusivity" /> 86 </file>87 88 <file id="file8" name_suffix="_icemod" description="ice variables" > 89 <field field_ref="ice_pres" />89 <field field_ref="w_masstr" name="vovematr" /> 90 </file> 91 92 <file id="file6" name_suffix="_icemod" description="ice variables" > 90 93 <field field_ref="snowthic_cea" name="snd" long_name="surface_snow_thickness" /> 91 94 <field field_ref="icethic_cea" name="sit" long_name="sea_ice_thickness" /> 92 <field field_ref="iceprod_cea" name="sip" long_name="sea_ice_thickness" /> 93 <field field_ref="ist_ipa" /> 94 <field field_ref="ioceflxb" /> 95 <field field_ref="uice_ipa" /> 96 <field field_ref="vice_ipa" /> 97 <field field_ref="utau_ice" /> 98 <field field_ref="vtau_ice" /> 99 <field field_ref="qsr_io_cea" /> 100 <field field_ref="qns_io_cea" /> 101 <field field_ref="snowpre" /> 102 </file> 103 95 <field field_ref="icevolu" name="iicevolu" /> 96 <field field_ref="snowvol" name="isnowvol" /> 97 <field field_ref="iceconc" name="iiceconc" /> 98 99 <field field_ref="icebopr" name="iicebopr" /> 100 <field field_ref="icedypr" name="iicedypr" /> 101 <field field_ref="ioceflxb" name="ioceflxb" /> 102 <field field_ref="uice_ipa" name="iicevelu" /> 103 <field field_ref="vice_ipa" name="iicevelv" /> 104 <field field_ref="isst" name="isstempe" /> 105 <field field_ref="isss" name="isssalin" /> 106 <field field_ref="micesalt" name="iicesali" /> 107 <field field_ref="miceage" name="iiceages" /> 108 <field field_ref="icelapr" name="iicelapr" /> 109 <field field_ref="icesipr" name="iicesipr" /> 110 <field field_ref="micet" name="iicetemp" /> 111 <field field_ref="icehc" name="iiceheco" /> 112 <field field_ref="isnowhc" name="isnoheco" /> 113 <field field_ref="icest" name="iicesurt" /> 114 <field field_ref="sfxbri" name="iicefsbr" /> 115 <field field_ref="sfxthd" name="iicefseq" /> 116 <field field_ref="ibrinv" name="ibrinvol" /> 117 <field field_ref="icecolf" name="iicecolf" /> 118 <field field_ref="icestr" name="iicestre" /> 119 <field field_ref="icevel" name="iicevelo" /> 120 <field field_ref="isume" name="iicesume" /> 121 <field field_ref="ibome" name="iicebome" /> 122 <field field_ref="idive" name="iicedive" /> 123 <field field_ref="ishear" name="iiceshea" /> 124 <field field_ref="icerepr" name="iicerepr" /> 125 <field field_ref="sfxmec" name="iicefsrp" /> 126 <field field_ref="sfxres" name="iicefsre" /> 127 <field field_ref="icetrp" name="iicevtrp" /> 128 </file> 129 130 131 <file id="file7" name_suffix="_scalar" description="scalar variables" > 132 <field field_ref="voltot" name="scvoltot" /> 133 <field field_ref="sshtot" name="scsshtot" /> 134 <field field_ref="sshsteric" name="scsshste" /> 135 <field field_ref="sshthster" name="scsshtst" /> 136 <field field_ref="masstot" name="scmastot" /> 137 <field field_ref="temptot" name="sctemtot" /> 138 <field field_ref="saltot" name="scsaltot" /> 139 140 <!-- available with ln_diahsb --> 141 <field field_ref="bgtemper" name="bgtemper" /> 142 <field field_ref="bgsaline" name="bgsaline" /> 143 <field field_ref="bgheatco" name="bgheatco" /> 144 <field field_ref="bgsaltco" name="bgsaltco" /> 145 <field field_ref="bgvolssh" name="bgvolssh" /> 146 <field field_ref="bgvoltot" name="bgvoltot" /> 147 <field field_ref="bgsshtot" name="bgsshtot" /> 148 <field field_ref="bgfrcvol" name="bgfrcvol" /> 149 <field field_ref="bgfrctem" name="bgfrctem" /> 150 <field field_ref="bgfrcsal" name="bgfrcsal" /> 151 <!-- available with ln_limdiahsb --> 152 <field field_ref="ibgvoltot" name="ibgvoltot" /> 153 <field field_ref="sbgvoltot" name="sbgvoltot" /> 154 <field field_ref="ibgarea" name="ibgarea" /> 155 <field field_ref="ibgsaline" name="ibgsaline" /> 156 <field field_ref="ibgtemper" name="ibgtemper" /> 157 <field field_ref="ibgheatco" name="ibgheatco" /> 158 <field field_ref="ibgsaltco" name="ibgsaltco" /> 159 <field field_ref="sbgheatco" name="sbgheatco" /> 160 <field field_ref="ibgfrcsfx" name="ibgfrcsfx" /> 161 <field field_ref="ibgfrcemp" name="ibgfrcemp" /> 162 <field field_ref="ibgsfx" name="ibgsfx" /> 163 <field field_ref="ibgemp" name="ibgemp" /> 164 <field field_ref="ibgsfxbri" name="ibgsfxbri" /> 165 <field field_ref="ibgsfxthd" name="ibgsfxthd" /> 166 <field field_ref="ibgsfxres" name="ibgsfxres" /> 167 <field field_ref="ibgsfxmec" name="ibgsfxmec" /> 168 <field field_ref="ibggrpme" name="ibggrpme" /> 169 170 </file> 171 172 <!-- 173 <file id="file8" name_suffix="_Tides" description="tidal harmonics" > 174 <field field_ref="M2x" name="M2x" long_name="M2 Elevation harmonic real part" /> 175 <field field_ref="M2y" name="M2y" long_name="M2 Elevation harmonic imaginary part" /> 176 <field field_ref="M2x_u" name="M2x_u" long_name="M2 current barotrope along i-axis harmonic real part " /> 177 <field field_ref="M2y_u" name="M2y_u" long_name="M2 current barotrope along i-axis harmonic imaginary part " /> 178 <field field_ref="M2x_v" name="M2x_v" long_name="M2 current barotrope along j-axis harmonic real part " /> 179 <field field_ref="M2y_v" name="M2y_v" long_name="M2 current barotrope along j-axis harmonic imaginary part " /> 180 </file> 181 --> 182 104 183 </file_group> 105 184 185 <file_group id="3d" output_freq="3d" output_level="10" enabled=".TRUE."/> <!-- 3d files --> 186 <file_group id="5d" output_freq="5d" output_level="10" enabled=".TRUE."/> <!-- 5d files --> 187 106 188 <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files --> 189 190 107 191 <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files --> 108 192 <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files --> … … 154 238 We must have buffer_size > jpi*jpj*jpk*8 (with jpi and jpj the subdomain size) 155 239 --> 156 <variable id="buffer_size" type="integer"> 25000000</variable>240 <variable id="buffer_size" type="integer">5000000</variable> 157 241 <variable id="buffer_server_factor_size" type="integer">2</variable> 158 242 <variable id="info_level" type="integer">0</variable> 159 <variable id="using_server" type="boolean"> false</variable>243 <variable id="using_server" type="boolean">true</variable> 160 244 <variable id="using_oasis" type="boolean">false</variable> 161 245 <variable id="oasis_codes_id" type="string" >oceanx</variable> -
branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/ORCA2_LIM3/cpp_ORCA2_LIM3.fcm
r4099 r4161 1 bld::tool::fppkeys key_trabbl key_ orca_r2 key_lim3 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_zdftke key_zdfddm key_zdftmx key_iomput key_mpp_mpi1 bld::tool::fppkeys key_trabbl key_lim3 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_zdftke key_zdfddm key_zdftmx key_iomput key_mpp_mpi -
branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/1_namelist_ref
r4153 r4161 281 281 ln_taudif = .false. ! HF tau contribution: use "mean of stress module - module of the mean stress" data 282 282 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow) 283 rn_vfac = 0. ! multiplicative factor for ocean/ice velocity 284 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.) 285 ! in the calculation of the wind stress (0.=absolute winds or 1.=relative winds) 283 286 / 284 287 !----------------------------------------------------------------------- … … 350 353 rn_si0 = 0.35 ! RGB & 2 bands: shortess depth of extinction 351 354 rn_si1 = 23.0 ! 2 bands: longest depth of extinction 355 ln_qsr_ice = .true. ! light penetration for ice-model LIM3 352 356 / 353 357 !----------------------------------------------------------------------- -
branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/field_def.xml
r4153 r4161 15 15 <!-- T grid --> 16 16 17 <field_group id="grid_T" grid_ref="grid_T_2D" >17 <field_group id="grid_T" grid_ref="grid_T_2D" > 18 18 <field id="toce" long_name="temperature" unit="degC" grid_ref="grid_T_3D"/> 19 19 <field id="soce" long_name="salinity" unit="psu" grid_ref="grid_T_3D"/> … … 50 50 </field_group> 51 51 52 <field_group id="Tides_T" grid_ref="grid_T_2D" operation="once" > 53 <!-- tidal composante --> 54 <field id="M2x" long_name="M2 Elevation harmonic real part " unit="m" /> 55 <field id="M2y" long_name="M2 Elevation harmonic imaginary part" unit="m" /> 56 <field id="S2x" long_name="S2 Elevation harmonic real part " unit="m" /> 57 <field id="S2y" long_name="S2 Elevation harmonic imaginary part" unit="m" /> 58 <field id="N2x" long_name="N2 Elevation harmonic real part " unit="m" /> 59 <field id="N2y" long_name="N2 Elevation harmonic imaginary part" unit="m" /> 60 <field id="K1x" long_name="K1 Elevation harmonic real part " unit="m" /> 61 <field id="K1y" long_name="K1 Elevation harmonic imaginary part" unit="m" /> 62 <field id="O1x" long_name="O1 Elevation harmonic real part " unit="m" /> 63 <field id="O1y" long_name="O1 Elevation harmonic imaginary part" unit="m" /> 64 <field id="Q1x" long_name="Q1 Elevation harmonic real part " unit="m" /> 65 <field id="Q1y" long_name="Q1 Elevation harmonic imaginary part" unit="m" /> 66 <field id="M4x" long_name="M4 Elevation harmonic real part " unit="m" /> 67 <field id="M4y" long_name="M4 Elevation harmonic imaginary part" unit="m" /> 68 <field id="K2x" long_name="K2 Elevation harmonic real part " unit="m" /> 69 <field id="K2y" long_name="K2 Elevation harmonic imaginary part" unit="m" /> 70 <field id="P1x" long_name="P1 Elevation harmonic real part " unit="m" /> 71 <field id="P1y" long_name="P1 Elevation harmonic imaginary part" unit="m" /> 72 <field id="Mfx" long_name="Mf Elevation harmonic real part " unit="m" /> 73 <field id="Mfy" long_name="Mf Elevation harmonic imaginary part" unit="m" /> 74 <field id="Mmx" long_name="Mm Elevation harmonic real part " unit="m" /> 75 <field id="Mmy" long_name="Mm Elevation harmonic imaginary part" unit="m" /> 76 </field_group> 77 78 <field_group id="Tides_U" grid_ref="grid_U_2D" operation="once" > 79 <field id="M2x_u" long_name="M2 current barotrope along i-axis harmonic real part " unit="m/s" /> 80 <field id="M2y_u" long_name="M2 current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 81 <field id="S2x_u" long_name="S2 current barotrope along i-axis harmonic real part " unit="m/s" /> 82 <field id="S2y_u" long_name="S2 current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 83 <field id="N2x_u" long_name="N2 current barotrope along i-axis harmonic real part " unit="m/s" /> 84 <field id="N2y_u" long_name="N2 current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 85 <field id="K1x_u" long_name="K1 current barotrope along i-axis harmonic real part " unit="m/s" /> 86 <field id="K1y_u" long_name="K1 current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 87 <field id="O1x_u" long_name="O1 current barotrope along i-axis harmonic real part " unit="m/s" /> 88 <field id="O1y_u" long_name="O1 current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 89 <field id="Q1x_u" long_name="Q1 current barotrope along i-axis harmonic real part " unit="m/s" /> 90 <field id="Q1y_u" long_name="Q1 current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 91 <field id="M4x_u" long_name="M4 current barotrope along i-axis harmonic real part " unit="m/s" /> 92 <field id="M4y_u" long_name="M4 current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 93 <field id="K2x_u" long_name="K2 current barotrope along i-axis harmonic real part " unit="m/s" /> 94 <field id="K2y_u" long_name="K2 current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 95 <field id="P1x_u" long_name="P1 current barotrope along i-axis harmonic real part " unit="m/s" /> 96 <field id="P1y_u" long_name="P1 current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 97 <field id="Mfx_u" long_name="Mf current barotrope along i-axis harmonic real part " unit="m/s" /> 98 <field id="Mfy_u" long_name="Mf current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 99 <field id="Mmx_u" long_name="Mm current barotrope along i-axis harmonic real part " unit="m/s" /> 100 <field id="Mmy_u" long_name="Mm current barotrope along i-axis harmonic imaginary part " unit="m/s" /> 101 </field_group> 102 103 <field_group id="Tides_V" grid_ref="grid_V_2D" operation="once" > 104 <field id="M2x_v" long_name="M2 current barotrope along j-axis harmonic real part " unit="m/s" /> 105 <field id="M2y_v" long_name="M2 current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 106 <field id="S2x_v" long_name="S2 current barotrope along j-axis harmonic real part " unit="m/s" /> 107 <field id="S2y_v" long_name="S2 current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 108 <field id="N2x_v" long_name="N2 current barotrope along j-axis harmonic real part " unit="m/s" /> 109 <field id="N2y_v" long_name="N2 current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 110 <field id="K1x_v" long_name="K1 current barotrope along j-axis harmonic real part " unit="m/s" /> 111 <field id="K1y_v" long_name="K1 current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 112 <field id="O1x_v" long_name="O1 current barotrope along j-axis harmonic real part " unit="m/s" /> 113 <field id="O1y_v" long_name="O1 current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 114 <field id="Q1x_v" long_name="Q1 current barotrope along j-axis harmonic real part " unit="m/s" /> 115 <field id="Q1y_v" long_name="Q1 current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 116 <field id="M4x_v" long_name="M4 current barotrope along j-axis harmonic real part " unit="m/s" /> 117 <field id="M4y_v" long_name="M4 current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 118 <field id="K2x_v" long_name="K2 current barotrope along j-axis harmonic real part " unit="m/s" /> 119 <field id="K2y_v" long_name="K2 current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 120 <field id="P1x_v" long_name="P1 current barotrope along j-axis harmonic real part " unit="m/s" /> 121 <field id="P1y_v" long_name="P1 current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 122 <field id="Mfx_v" long_name="Mf current barotrope along j-axis harmonic real part " unit="m/s" /> 123 <field id="Mfy_v" long_name="Mf current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 124 <field id="Mmx_v" long_name="Mm current barotrope along j-axis harmonic real part " unit="m/s" /> 125 <field id="Mmy_v" long_name="Mm current barotrope along j-axis harmonic imaginary part " unit="m/s" /> 126 </field_group> 127 128 52 129 <!-- SBC --> 53 130 … … 59 136 <field id="snowpre" long_name="Snow precipitation" unit="kg/m2/s" /> 60 137 <field id="runoffs" long_name="River Runoffs" unit="Kg/m2/s" /> 138 <field id="precip" long_name="Total precipitation" unit="kg/m2/s" /> 139 61 140 62 141 <field id="qt" long_name="Net Downward Heat Flux" unit="W/m2" /> … … 79 158 <field id="qhc_oce" long_name="Downward Heat Content of E-P over open ocean" unit="W/m2" /> 80 159 <field id="taum_oce" long_name="wind stress module over open ocean" unit="N/m2" /> 81 82 <field id="ice_cover" long_name="Ice fraction" unit="1" />83 84 <field id="ioceflxb" long_name="Oceanic flux at the ice base" unit="W/m2" />85 <field id="qsr_ai_cea" long_name="Air-Ice downward solar heat flux (cell average)" unit="W/m2" />86 <field id="qns_ai_cea" long_name="Air-Ice downward non-solar heat flux (cell average)" unit="W/m2" />87 <field id="qla_ai_cea" long_name="Air-Ice downward Latent heat flux (cell average)" unit="W/m2" />88 89 <field id="qsr_io_cea" long_name="Ice-Oce downward solar heat flux (cell average)" unit="W/m2" />90 <field id="qns_io_cea" long_name="Ice-Oce downward non-solar heat flux (cell average)" unit="W/m2" />91 92 <field id="snowthic_cea" long_name="Snow thickness (cell average)" unit="m" />93 <field id="icethic_cea" long_name="Ice thickness (cell average)" unit="m" />94 <field id="iceprod_cea" long_name="Ice production (cell average)" unit="m/s" />95 96 <field id="ice_pres" long_name="Ice presence" unit="-" />97 <field id="ist_cea" long_name="Ice surface temperature (cell average)" unit="degC" />98 <field id="ist_ipa" long_name="Ice surface temperature (ice presence average)" unit="degC" />99 <field id="uice_ipa" long_name="Ice velocity along i-axis at I-point (ice presence average)" unit="m/s" />100 <field id="vice_ipa" long_name="Ice velocity along j-axis at I-point (ice presence average)" unit="m/s" />101 102 <field id="utau_ice" long_name="Wind stress along i-axis over the ice at i-point" unit="N/m2" />103 <field id="vtau_ice" long_name="Wind stress along j-axis over the ice at i-point" unit="N/m2" />104 105 <field id="u_imasstr" long_name="Sea-ice mass transport along i-axis" unit="kg/s" />106 <field id="v_imasstr" long_name="Sea-ice mass transport along j-axis" unit="kg/s" />107 <field id="emp_x_sst" long_name="Concentration/Dilution term on SST" unit="kgC/m2/s" />108 <field id="emp_x_sss" long_name="Concentration/Dilution term on SSS" unit="kgPSU/m2/s" />109 160 110 161 <!-- available key_coupled --> … … 132 183 <field id="sntoice_cea" long_name="Snow-Ice Formation Rate (cell average)" unit="kg/m2/s" /> 133 184 <field id="ticemel_cea" long_name="Rate of Melt at Upper Surface of Sea Ice (cell average)" unit="kg/m2/s" /> 185 186 <!-- ice fields --> 187 188 <field id="ice_cover" long_name="Ice fraction" unit="1" /> 189 190 <field id="ioceflxb" long_name="Oceanic flux at the ice base" unit="W/m2" /> 191 <field id="qsr_ai_cea" long_name="Air-Ice downward solar heat flux (cell average)" unit="W/m2" /> 192 <field id="qns_ai_cea" long_name="Air-Ice downward non-solar heat flux (cell average)" unit="W/m2" /> 193 <field id="qla_ai_cea" long_name="Air-Ice downward Latent heat flux (cell average)" unit="W/m2" /> 194 195 <field id="qsr_io_cea" long_name="Ice-Oce downward solar heat flux (cell average)" unit="W/m2" /> 196 <field id="qns_io_cea" long_name="Ice-Oce downward non-solar heat flux (cell average)" unit="W/m2" /> 197 198 <field id="snowthic_cea" long_name="Snow thickness (cell average)" unit="m" /> 199 <field id="icethic_cea" long_name="Ice thickness (cell average)" unit="m" /> 200 <field id="iceprod_cea" long_name="Ice production (cell average)" unit="m/s" /> 201 202 <field id="ice_pres" long_name="Ice presence" unit="-" /> 203 <field id="ist_cea" long_name="Ice surface temperature (cell average)" unit="degC" /> 204 <field id="ist_ipa" long_name="Ice surface temperature (ice presence average)" unit="degC" /> 205 <field id="uice_ipa" long_name="Ice velocity along i-axis at I-point (ice presence average)" unit="m/s" /> 206 <field id="vice_ipa" long_name="Ice velocity along j-axis at I-point (ice presence average)" unit="m/s" /> 207 208 <field id="utau_ice" long_name="Wind stress along i-axis over the ice at i-point" unit="N/m2" /> 209 <field id="vtau_ice" long_name="Wind stress along j-axis over the ice at i-point" unit="N/m2" /> 210 211 <field id="u_imasstr" long_name="Sea-ice mass transport along i-axis" unit="kg/s" /> 212 <field id="v_imasstr" long_name="Sea-ice mass transport along j-axis" unit="kg/s" /> 213 <field id="emp_x_sst" long_name="Concentration/Dilution term on SST" unit="kgC/m2/s" /> 214 <field id="emp_x_sss" long_name="Concentration/Dilution term on SSS" unit="kgPSU/m2/s" /> 215 216 217 <field id="iceconc" long_name="ice concentration" unit="%" /> 218 <field id="icebopr" long_name="daily bottom thermo ice prod." unit="km3/day" /> 219 <field id="icedypr" long_name="daily dynamic ice prod." unit="km3/day" /> 220 <field id="ioceflxb" long_name="Oceanic flux at the ice base" unit="W/m2" /> 221 <field id="uice_ipa" long_name="Ice velocity along i-axis at I-point (ice presence average)" unit="m/s" /> 222 <field id="vice_ipa" long_name="Ice velocity along j-axis at I-point (ice presence average)" unit="m/s" /> 223 <field id="isst" long_name="sea surface temperature" unit="degC" /> 224 <field id="isss" long_name="sea surface salinity" unit="psu" /> 225 <field id="qt_oce" long_name="total flux at ocean surface" unit="W/m2" /> 226 <field id="qsr_oce" long_name="solar heat flux at ocean surface" unit="W/m2" /> 227 <field id="qns_oce" long_name="non-solar heat flux at ocean surface" unit="W/m2" /> 228 <field id="hfbri" long_name="heat flux due to brine release" unit="W/m2" /> 229 <field id="utau_ice" long_name="Wind stress along i-axis over the ice at i-point" unit="N/m2" /> 230 <field id="vtau_ice" long_name="Wind stress along j-axis over the ice at i-point" unit="N/m2" /> 231 <field id="qsr_io" long_name="Ice-Oce downward solar heat flux" unit="W/m2" /> 232 <field id="qns_io" long_name="Ice-Oce downward non-solar heat flux" unit="W/m2" /> 233 <field id="micesalt" long_name="Mean ice salinity" unit="psu" /> 234 <field id="miceage" long_name="Mean ice age" unit="years" /> 235 <field id="icelapr" long_name="daily lateral thermo ice prod." unit="km3/day" /> 236 <field id="icesipr" long_name="daily snowice ice prod." unit="km3/day" /> 237 <field id="micet" long_name="Mean ice temperature" unit="degC" /> 238 <field id="icehc" long_name="ice total heat content" unit="10^9 J" /> 239 <field id="isnowhc" long_name="snow total heat content" unit="10^9J" /> 240 <field id="icest" long_name="ice surface temperature" unit="degC" /> 241 <field id="sfxbri" long_name="brine salt flux" unit="psu*kg/m2/day" /> 242 <field id="sfxthd" long_name="equivalent FW salt flux" unit="psu*kg/m2/day" /> 243 <field id="ibrinv" long_name="brine volume" unit="%" /> 244 <field id="icecolf" long_name="frazil ice collection thickness" unit="m" /> 245 <field id="icestr" long_name="ice strength" unit="N/m" /> 246 <field id="icevel" long_name="ice velocity" unit="m/s" /> 247 <field id="isume" long_name="surface melt" unit="km3/day" /> 248 <field id="ibome" long_name="bottom melt" unit="km3/day" /> 249 <field id="idive" long_name="divergence" unit="10-8s-1" /> 250 <field id="ishear" long_name="shear" unit="10-8s-1" /> 251 <field id="icerepr" long_name="daily resultant ice prod./melting from limupdate" unit="km3/day" /> 252 <field id="icevolu" long_name="ice volume" unit="km3" /> 253 <field id="snowvol" long_name="snow volume" unit="km3" /> 254 <field id="sfxmec" long_name="salt flux from ridging rafting" unit="psu*kg/m2/day" /> 255 <field id="sfxres" long_name="salt flux from lipupdate (resultant)" unit="psu*kg/m2/day" /> 256 <field id="icetrp" long_name="ice volume transport" unit="km3/day" /> 257 134 258 135 259 </field_group> … … 214 338 <field id="saltot" long_name="global mean salinity" unit="psu" /> 215 339 <field id="fram_trans" long_name="Sea Ice Mass Transport Through Fram Strait" unit="kg/s" /> 340 <!-- available with ln_diahsb --> 341 <field id="bgtemper" long_name="global mean temperature" unit="degC" /> 342 <field id="bgsaline" long_name="global mean salinity" unit="psu" /> 343 <field id="bgheatco" long_name="global mean heat content" unit="10^9J" /> 344 <field id="bgsaltco" long_name="global mean salt content" unit="psu*m3" /> 345 <field id="bgvolssh" long_name="global mean ssh volume" unit="km3" /> 346 <field id="bgvoltot" long_name="global mean volume" unit="km3" /> 347 <field id="bgsshtot" long_name="global mean ssh" unit="m" /> 348 <field id="bgfrcvol" long_name="global mean volume from forcing" unit="km3" /> 349 <field id="bgfrctem" long_name="global mean heat content from forcing" unit="10^9J" /> 350 <field id="bgfrcsal" long_name="global mean salt content from forcing" unit="psu*km3" /> 351 </field_group> 352 353 <field_group id="SBC_scalar" domain_ref="1point" > 354 <!-- available with ln_limdiahsb --> 355 <field id="ibgvoltot" long_name="global mean ice volume" unit="km3" /> 356 <field id="sbgvoltot" long_name="global mean snow volume" unit="km3" /> 357 <field id="ibgarea" long_name="global mean ice area" unit="km2" /> 358 <field id="ibgsaline" long_name="global mean ice salinity" unit="psu" /> 359 <field id="ibgtemper" long_name="global mean ice temperature" unit="degC" /> 360 <field id="ibgheatco" long_name="global mean ice heat content" unit="10^9J" /> 361 <field id="ibgsaltco" long_name="global mean ice salt content" unit="psu*km3" /> 362 <field id="sbgheatco" long_name="global mean snow heat content" unit="10^9J" /> 363 <field id="ibgfrcsfx" long_name="global mean salt content from sfx" unit="psu*km3" /> 364 <field id="ibgfrcemp" long_name="global mean volume from emp" unit="km3" /> 365 <field id="ibgsfx" long_name="global mean emps" unit="psu*kg/m2/day" /> 366 <field id="ibgemp" long_name="global mean emp" unit="kg/m2/day" /> 367 <field id="ibgsfxbri" long_name="global mean ice sfx_bri" unit="psu*kg/m2/day" /> 368 <field id="ibgsfxthd" long_name="global mean ice sfx_thd" unit="psu*kg/m2/day" /> 369 <field id="ibgsfxres" long_name="global mean ice sfx_res" unit="psu*kg/m2/day" /> 370 <field id="ibgsfxmec" long_name="global mean ice fsalt_rpo" unit="psu*kg/m2/day" /> 371 <field id="ibggrpme" long_name="global mean ice growth+melt volume" unit="km3" /> 216 372 </field_group> 217 373 -
branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ice_lim2_ref
r4148 r4161 54 54 telast = 9600 ! timescale for EVP elastic waves 55 55 alphaevp = 1.0 ! coefficient for the solution of EVP int. stresses 56 hminrhg = 0.05 ! ice thickness (m) below which ice velocity equal ocean velocity 56 57 / 57 58 !----------------------------------------------------------------------- -
branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r4147 r4161 14 14 &namicerun ! Share parameters for dynamics/advection/thermo 15 15 !----------------------------------------------------------------------- 16 cn_icerst_in = "restart_ice _in" ! suffix of ice restart name (input)16 cn_icerst_in = "restart_ice" ! suffix of ice restart name (input) 17 17 cn_icerst_out = "restart_ice" ! suffix of ice restart name (output) 18 18 ln_limdyn = .true. ! ice dynamics (T) or thermodynamics only (F) 19 acrit = 1.0e-02 , 1.0e-02 ! minimum fraction for leads in the Northern (Southern) Hemisphere 20 hsndif = 0.0 ! computation of temperature in snow (=0.0) or not 21 hicdif = 0.0 ! computation of temperature in ice (=0.0) or not (=9999.0) 19 amax = 0.999 ! maximum ice concentration 22 20 cai = 1.40e-3 ! atmospheric drag over sea ice 23 21 cao = 1.00e-3 ! atmospheric drag over ocean 24 22 ln_nicep = .false. ! Ice points output for debug (yes or no) 23 ln_limdiahsb = .false. ! check the heat and salt budgets (T) or not (F) 24 ln_limdiaout = .false. ! output the heat and salt budgets (T) or not (F) 25 25 / 26 26 !----------------------------------------------------------------------- … … 30 30 hninn = 0.3 ! initial snow thickness in the north 31 31 hnins = 0.1 ! " " south 32 hginn_u = 3.50 ! initial undeformed ice thickness in the north 33 hgins_u = 1.0 ! " " " south 34 aginn_u = 0.95 ! initial undeformed ice concentration in the north 35 agins_u = 0.9 ! " " " south 36 hginn_d = 0.0 ! initial deformed ice thickness in the north 37 hgins_d = 0.0 ! 38 aginn_d = 0.00 ! initial deformed ice concentration in the north 39 agins_d = 0.00 ! " " " south 32 hginn = 3.50 ! initial undeformed ice thickness in the north 33 hgins = 1.0 ! " " " south 34 aginn = 0.95 ! initial undeformed ice concentration in the north 35 agins = 0.9 ! " " " south 40 36 sinn = 6.301 ! initial salinity in the north 41 37 sins = 6.301 ! " " south … … 62 58 telast =9600.0 ! timescale for elastic waves, SB, 720.0 63 59 alphaevp = 1.0 ! coefficient for the solution of internal ice stresses 60 hminrhg = 0.05 ! ice thickness (m) below which ice velocity equal ocean velocity 64 61 / 65 62 !----------------------------------------------------------------------- … … 80 77 hicmin = 0.2 ! ice thickness corr. to max. energy stored in brine pocket 81 78 hiclim = 0.10 ! minimum ice thickness 82 amax = 0.999 ! maximum lead fraction83 79 sbeta = 1. ! numerical caracteritic of the scheme for diffusion in ice 84 80 ! Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) … … 145 141 ninfo = 1 ! frequency of ouputs on file ice_evolu in case of averaging 146 142 / 143 !!----------------------------------------------------------------------- 144 !&namicehsb ! Heat and salt budgets 145 !!----------------------------------------------------------------------- 146 ! 147 !/ 147 148 !----------------------------------------------------------------------- 148 149 &namiceout ! parameters for outputs 149 150 !----------------------------------------------------------------------- 150 noumef = 37! number of fields151 noumef = 43 ! number of fields 151 152 add_diag_swi= 1 ! 1 -> diagnose distribution in thickness space 152 153 ! 0 -> only simple diagnostics … … 157 158 field_2 = 'Ice thickness ', 'iicethic', 'm ', 1 , 1.0 , 0.0 158 159 field_3 = 'Snow thickness ', 'isnowthi', 'm ', 1 , 1.0 , 0.0 159 field_4 = 'Daily bottom thermo ice production ', 'iicebopr', ' cm/day ', 1 , 100., 0.0160 field_5 = 'Daily dynamic ice production ', 'iicedypr', ' cm/day ', 1 , 100., 0.0160 field_4 = 'Daily bottom thermo ice production ', 'iicebopr', 'km3/day ', 1 , 1.0e-9 , 0.0 161 field_5 = 'Daily dynamic ice production ', 'iicedypr', 'km3/day ', 1 , 1.0e-9 , 0.0 161 162 field_6 = 'Oceanic flux at the ice base ', 'ioceflxb', 'w/m2 ', 1 , 1.0 , 0.0 162 163 field_7 = 'Ice velocity u ', 'iicevelu', 'm/s ', 1 , 1.0 , 0.0 … … 172 173 field_17 = 'Solar flux at ice/ocean surface ', 'iicesflx', 'w/m2 ', 1 , 1.0 , 0.0 173 174 field_18 = 'Non-solar flux at ice/ocean surface', 'iicenflx', 'w/m2 ', 1 , 1.0 , 0.0 174 field_19 = 'Snow precipitation ', 'isnowpre', 'kg/ day', 1 , 1.0 , 0.0175 field_19 = 'Snow precipitation ', 'isnowpre', 'kg/m2/d ', 1 , 1.0 , 0.0 175 176 field_20 = 'Mean ice salinity ', 'iicesali', 'psu ', 1 , 1.0 , 0.0 176 177 field_21 = 'Mean ice age ', 'iiceages', 'years ', 1 , 0.002739, 0.0 177 field_22 = 'Daily lateral thermo ice prod. ', 'iicelapr', ' cm/day ', 1 ,100., 0.0178 field_23 = 'Daily snowice ice production ', 'iicesipr', ' cm/day ', 1 ,100., 0.0178 field_22 = 'Daily lateral thermo ice prod. ', 'iicelapr', 'km3/day ', 1 ,1.0e-9 , 0.0 179 field_23 = 'Daily snowice ice production ', 'iicesipr', 'km3/day ', 1 ,1.0e-9 , 0.0 179 180 field_24 = 'Mean ice temperature ', 'iicetemp', 'C ', 1 , 1.0 , -273.15 180 181 field_25 = 'Ice total heat content ', 'iiceheco', '10^9 J ', 1 , 1.0 , 0.0 181 182 field_26 = 'Ice surface temperature ', 'iicesurt', 'C ', 1 , 1.0 , -273.15 182 183 field_27 = 'Snow temperature ', 'isnotem2', 'C ', 1 , 1.0 , -273.15 183 field_28 = 'Fsbri - brine salt flux ', 'iic fsbri', 'kg/m2/s', 1 , 1.0 , 0.0184 field_29 = 'Fseqv - equivalent FW salt flux ', 'iic fseqv', 'kg/m2/s', 1 , 1.0 , 0.0184 field_28 = 'Fsbri - brine salt flux ', 'iicefsbr', 'kg/m2/d ', 1 , 1.0 , 0.0 185 field_29 = 'Fseqv - equivalent FW salt flux ', 'iicefseq', 'kg/m2/d ', 1 , 1.0 , 0.0 185 186 field_30 = 'Brine volume ', 'ibrinvol', '% ', 1 , 100.0 , 0.0 186 187 field_31 = 'Frazil ice collection thickness ', 'iicecolf', 'm ', 1 , 1.0 , 0.0 187 188 field_32 = 'Ice strength ', 'iicestre', 'N/m ', 1 , 0.001 , 0.0 188 189 field_33 = 'Ice velocity ', 'iicevelo', 'm/s ', 1 , 1.0 , 0.0 189 field_34 = 'Surface melt ', 'iicesume', ' cm/day ', 1 ,100., 0.0190 field_35 = 'Bottom melt ', 'iicebome', ' cm/day ', 1 ,100., 0.0190 field_34 = 'Surface melt ', 'iicesume', 'km3/day ', 1 ,1.0e-9 , 0.0 191 field_35 = 'Bottom melt ', 'iicebome', 'km3/day ', 1 ,1.0e-9 , 0.0 191 192 field_36 = 'Divergence ', 'iicedive', '10-8s-1 ', 1 , 1.0e8 , 0.0 192 193 field_37 = 'Shear ', 'iiceshea', '10-8s-1 ', 1 , 1.0e8 , 0.0 193 / 194 field_38 = 'Daily resultant ice prod/melt ', 'iicerepr', 'km3/day ', 1 , 1.0e-9 , 0.0 195 field_39 = 'Ice volume ', 'iicevolu', 'km3 ', 1 , 1.0e-9 , 0.0 196 field_40 = 'Snow volume ', 'isnowvol', 'km3 ', 1 , 1.0e-9 , 0.0 197 field_41 = 'Fsrpo - salt flux from ridg/raft ', 'iicefsrp', 'kg/m2/d ', 1 , 1.0 , 0.0 198 field_42 = 'Fsres - salt flux from limupdate ', 'iicefsre', 'kg/m2/d ', 1 , 1.0 , 0.0 199 field_43 = 'Ice volume transport ', 'iicevtrp', 'km3/day ', 1 ,1.0e-9 , 0.0 200 / 201 -
branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ref
r4152 r4161 283 283 ln_taudif = .false. ! HF tau contribution: use "mean of stress module - module of the mean stress" data 284 284 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow) 285 rn_vfac = 0. ! multiplicative factor for ocean/ice velocity 286 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.) 287 ! in the calculation of the wind stress (0.=absolute winds or 1.=relative winds) 285 288 / 286 289 !----------------------------------------------------------------------- … … 352 355 rn_si0 = 0.35 ! RGB & 2 bands: shortess depth of extinction 353 356 rn_si1 = 23.0 ! 2 bands: longest depth of extinction 357 ln_qsr_ice = .true. ! light penetration for ice-model LIM3 354 358 / 355 359 !----------------------------------------------------------------------- -
branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/cfg.txt
r4159 r4161 1 ORCA2_LIM3 OPA_SRC LIM_SRC_3 1 GYRE OPA_SRC 2 2 GYRE_BFM OPA_SRC TOP_SRC 3 GYRE OPA_SRC3 GYRE_PISCES OPA_SRC TOP_SRC 4 4 AMM12 OPA_SRC 5 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 6 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 7 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 8 ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 5 9 ORCA2_SAS_LIM OPA_SRC SAS_SRC LIM_SRC_2 NST_SRC 6 ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC7 GYRE_PISCES OPA_SRC TOP_SRC8 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC9 10 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90
r4147 r4161 51 51 REAL(wp), PUBLIC :: ahi0 !: sea-ice hor. eddy diffusivity coeff. (m2/s) 52 52 REAL(wp), PUBLIC :: alphaevp !: coefficient for the solution of EVP int. stresses 53 REAL(wp), PUBLIC :: hminrhg = 0.001_wp !: clem : ice volume (a*h in m) below which ice velocity is set to ocean velocity 53 54 54 55 REAL(wp), PUBLIC :: usecc2 !: = 1.0 / ( ecc * ecc ) -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90
r4147 r4161 228 228 & dm, nbiter, nbitdr, om, resl, cw, angvg, pstar, & 229 229 & c_rhg, etamn, creepl, ecc, ahi0, & 230 & nevp, telast, alphaevp230 & nevp, telast, alphaevp, hminrhg 231 231 !!------------------------------------------------------------------- 232 232 … … 262 262 WRITE(numout,*) ' timescale for elastic waves telast = ', telast 263 263 WRITE(numout,*) ' coefficient for the solution of int. stresses alphaevp = ', alphaevp 264 WRITE(numout,*) ' min ice thickness for rheology calculations hminrhg = ', hminrhg 264 265 ENDIF 265 266 ! -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90
r3625 r4161 31 31 32 32 !!---------------------------------------------------------------------- 33 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)33 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 34 34 !! $Id$ 35 35 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r4147 r4161 188 188 REAL(wp), PUBLIC :: alphaevp !: coeficient of the internal stresses !SB 189 189 REAL(wp), PUBLIC :: unit_fac = 1.e+09_wp !: conversion factor for ice / snow enthalpy 190 REAL(wp), PUBLIC :: hminrhg = 0.001_wp !: clem : ice volume (a*h, in m) below which ice velocity is set to ocean velocity 190 191 191 192 ! !!** ice-salinity namelist (namicesal) ** … … 405 406 !!-------------------------------------------------------------------------- 406 407 !! Check if everything down here is necessary 408 LOGICAL , PUBLIC :: ln_limdiahsb !: flag for ice diag (T) or not (F) 409 LOGICAL , PUBLIC :: ln_limdiaout !: flag for ice diag (T) or not (F) 407 410 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_newice !: volume of ice formed in the leads 408 411 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dv_dt_thd !: thermodynamic growth rates … … 414 417 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_bot_me ! vertical bottom melt 415 418 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_sur_me ! vertical surface melt 419 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_res_pr ! production (growth+melt) due to limupdate 420 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vi ! transport of ice volume 416 421 INTEGER , PUBLIC :: jiindx, jjindx !: indexes of the debugging point 417 422 418 423 !!---------------------------------------------------------------------- 419 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2010)424 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 420 425 !! $Id$ 421 426 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 527 532 & izero (jpi,jpj,jpl) , diag_bot_gr(jpi,jpj) , diag_dyn_gr(jpi,jpj) , & 528 533 & fstroc (jpi,jpj,jpl) , diag_bot_me(jpi,jpj) , diag_sur_me(jpi,jpj) , & 529 & fhbricat (jpi,jpj,jpl) , v_newice (jpi,jpj), STAT=ierr(ii) )534 & fhbricat (jpi,jpj,jpl) , diag_res_pr(jpi,jpj) , diag_trp_vi(jpi,jpj) , v_newice(jpi,jpj) , STAT=ierr(ii) ) 530 535 531 536 ice_alloc = MAXVAL( ierr(:) ) -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
r4147 r4161 40 40 41 41 !!---------------------------------------------------------------------- 42 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)42 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 43 43 !! $Id$ 44 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 128 128 !! ** input : Namelist namicerun 129 129 !!------------------------------------------------------------------- 130 NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, a crit, hsndif, hicdif, cai, cao, ln_nicep130 NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb, ln_limdiaout 131 131 INTEGER :: ios ! Local integer output status for namelist read 132 132 !!------------------------------------------------------------------- … … 144 144 ln_nicep = .FALSE. 145 145 CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' ) 146 ENDIF 146 ENDIF 147 147 ! 148 148 IF(lwp) THEN ! control print … … 151 151 WRITE(numout,*) ' ~~~~~~' 152 152 WRITE(numout,*) ' switch for ice dynamics (1) or not (0) ln_limdyn = ', ln_limdyn 153 WRITE(numout,*) ' minimum fraction for leads in the NH (SH) acrit(1/2) = ', acrit(:) 154 WRITE(numout,*) ' computation of temp. in snow (=0) or not (=9999) hsndif = ', hsndif 155 WRITE(numout,*) ' computation of temp. in ice (=0) or not (=9999) hicdif = ', hicdif 153 WRITE(numout,*) ' maximum ice concentration = ', amax 156 154 WRITE(numout,*) ' atmospheric drag over sea ice = ', cai 157 155 WRITE(numout,*) ' atmospheric drag over ocean = ', cao 158 156 WRITE(numout,*) ' Several ice points in the ice or not in ocean.output = ', ln_nicep 157 WRITE(numout,*) ' Diagnose heat/salt budget or not ln_limdiahsb = ', ln_limdiahsb 158 WRITE(numout,*) ' Output heat/salt budget or not ln_limdiaout = ', ln_limdiaout 159 159 ENDIF 160 160 ! -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
r3764 r4161 15 15 !! lim_adv_y : advection of sea ice on y axis 16 16 !!---------------------------------------------------------------------- 17 USE dom_oce ! ocean domain 18 USE ice ! LIM-3 variables 19 USE dom_ice ! LIM-3 domain 20 USE lbclnk ! lateral boundary condition - MPP exchanges 21 USE in_out_manager ! I/O manager 22 USE prtctl ! Print control 23 USE lib_mpp ! MPP library 24 USE wrk_nemo ! work arrays 25 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 26 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 17 USE dom_oce ! ocean domain 18 USE dom_ice ! LIM-3 domain 19 USE ice ! LIM-3 variables 20 USE lbclnk ! lateral boundary condition - MPP exchanges 21 USE in_out_manager ! I/O manager 22 USE prtctl ! Print control 23 USE lib_mpp ! MPP library 24 USE wrk_nemo ! work arrays 25 USE lib_fortran ! to use key_nosignedzero 27 26 28 27 IMPLICIT NONE … … 39 38 # include "vectopt_loop_substitute.h90" 40 39 !!---------------------------------------------------------------------- 41 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)40 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 42 41 !! $Id$ 43 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r3625 r4161 30 30 31 31 !!---------------------------------------------------------------------- 32 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)32 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 33 33 !! $Id$ 34 34 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r4147 r4161 15 15 !! lim_dyn_init : initialization and namelist read 16 16 !!---------------------------------------------------------------------- 17 USE phycst ! physical constants 18 USE dom_oce ! ocean space and time domain 19 USE sbc_oce ! Surface boundary condition: ocean fields 20 USE sbc_ice ! Surface boundary condition: ice fields 21 USE ice ! LIM-3 variables 22 USE par_ice ! LIM-3 parameters 23 USE dom_ice ! LIM-3 domain 24 USE limrhg ! LIM-3 rheology 25 USE lbclnk ! lateral boundary conditions - MPP exchanges 26 USE lib_mpp ! MPP library 27 USE wrk_nemo ! work arrays 28 USE in_out_manager ! I/O manager 29 USE prtctl ! Print control 30 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 17 USE phycst ! physical constants 18 USE dom_oce ! ocean space and time domain 19 USE sbc_oce ! Surface boundary condition: ocean fields 20 USE sbc_ice ! Surface boundary condition: ice fields 21 USE ice ! LIM-3 variables 22 USE par_ice ! LIM-3 parameters 23 USE dom_ice ! LIM-3 domain 24 USE limrhg ! LIM-3 rheology 25 USE lbclnk ! lateral boundary conditions - MPP exchanges 26 USE lib_mpp ! MPP library 27 USE wrk_nemo ! work arrays 28 USE in_out_manager ! I/O manager 29 USE prtctl ! Print control 30 USE lib_fortran ! glob_sum 31 USE timing ! Timing 31 32 32 33 IMPLICIT NONE … … 38 39 # include "vectopt_loop_substitute.h90" 39 40 !!---------------------------------------------------------------------- 40 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)41 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 41 42 !! $Id$ 42 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 65 66 REAL(wp), POINTER, DIMENSION(:) :: zmsk ! i-averaged of tmask 66 67 REAL(wp), POINTER, DIMENSION(:,:) :: zu_io, zv_io ! ice-ocean velocity 67 !!--------------------------------------------------------------------- 68 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 69 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 70 !!--------------------------------------------------------------------- 71 72 IF( nn_timing == 1 ) CALL timing_start('limdyn') 68 73 69 74 CALL wrk_alloc( jpi, jpj, zu_io, zv_io ) 70 75 CALL wrk_alloc( jpj, zind, zmsk ) 76 77 ! ------------------------------- 78 !- check conservation (C Rousset) 79 IF (ln_limdiahsb) THEN 80 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 81 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 82 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 83 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 84 ENDIF 85 !- check conservation (C Rousset) 86 ! ------------------------------- 71 87 72 88 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) … … 208 224 ENDIF 209 225 ! 226 ! ------------------------------- 227 !- check conservation (C Rousset) 228 IF (ln_limdiahsb) THEN 229 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 230 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 231 232 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 233 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 234 235 zchk_vmin = glob_min(v_i) 236 zchk_amax = glob_max(SUM(a_i,dim=3)) 237 zchk_amin = glob_min(a_i) 238 239 IF(lwp) THEN 240 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limdyn) = ',(zchk_v_i * rday) 241 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * rday) 242 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limdyn) = ',(zchk_vmin * 1.e-3) 243 !IF ( zchk_amax > amax+1.e-10 ) WRITE(numout,*) 'violation a_i>amax (limdyn) = ',zchk_amax 244 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limdyn) = ',zchk_amin 245 ENDIF 246 ENDIF 247 !- check conservation (C Rousset) 248 ! ------------------------------- 249 210 250 CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 211 251 CALL wrk_dealloc( jpj, zind, zmsk ) 212 252 ! 253 IF( nn_timing == 1 ) CALL timing_stop('limdyn') 254 213 255 END SUBROUTINE lim_dyn 214 256 … … 230 272 & dm, nbiter, nbitdr, om, resl, cw, angvg, pstar, & 231 273 & c_rhg, etamn, creepl, ecc, ahi0, & 232 & nevp, telast, alphaevp 274 & nevp, telast, alphaevp, hminrhg 233 275 !!------------------------------------------------------------------- 234 276 … … 264 306 WRITE(numout,*) ' timescale for elastic waves telast = ', telast 265 307 WRITE(numout,*) ' coefficient for the solution of int. stresses alphaevp = ', alphaevp 308 WRITE(numout,*) ' min ice thickness for rheology calculations hminrhg = ', hminrhg 266 309 ENDIF 267 310 ! -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r3625 r4161 35 35 # include "vectopt_loop_substitute.h90" 36 36 !!---------------------------------------------------------------------- 37 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2010)37 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 38 38 !! $Id$ 39 39 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 137 137 END DO ! end of sub-time step loop 138 138 139 ! ----------------------- 140 !!! final step (clem) !!! 141 DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction 142 DO ji = 1 , fs_jpim1 ! vector opt. 143 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 144 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 145 END DO 146 END DO 147 ! 148 DO jj= 2, jpjm1 ! diffusive trend : divergence of the fluxes 149 DO ji = fs_2 , fs_jpim1 ! vector opt. 150 zdiv (ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj ) & 151 & + zflv(ji,jj) - zflv(ji ,jj-1) ) / ( e1t (ji,jj) * e2t (ji,jj) ) 152 ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) ) 153 END DO 154 END DO 155 CALL lbc_lnk( ptab, 'T', 1. ) ! lateral boundary condition 156 !!! final step (clem) !!! 157 ! ----------------------- 158 139 159 IF(ln_ctl) THEN 140 160 zrlx(:,:) = ptab(:,:) - ztab0(:,:) -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r4147 r4161 5 5 !!====================================================================== 6 6 !! History : 2.0 ! 2004-01 (C. Ethe, G. Madec) Original code 7 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 7 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 8 !! - ! 2012 (C. Rousset) add par_oce (for jp_sal)...bug? 8 9 !!---------------------------------------------------------------------- 9 10 #if defined key_lim3 … … 18 19 USE dom_oce ! ocean domain 19 20 USE sbc_oce ! Surface boundary condition: ocean fields 21 USE sbc_ice ! Surface boundary condition: ice fields 20 22 USE eosbn2 ! equation of state 21 23 USE ice ! sea-ice variables 22 24 USE par_ice ! ice parameters 25 USE par_oce ! ocean parameters 23 26 USE dom_ice ! sea-ice domain 24 27 USE in_out_manager ! I/O manager 25 28 USE lbclnk ! lateral boundary condition - MPP exchanges 26 29 USE lib_mpp ! MPP library 30 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 27 31 USE wrk_nemo ! work arrays 28 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)29 32 30 33 IMPLICIT NONE … … 49 52 50 53 !!---------------------------------------------------------------------- 51 !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)54 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) 52 55 !! $Id$ 53 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 54 !!---------------------------------------------------------------------- 56 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 57 !!---------------------------------------------------------------------- 58 55 59 CONTAINS 56 60 … … 61 65 !! ** Purpose : defined the sea-ice initial state 62 66 !! 63 !! ** Method : restart from a state defined in a binary file 64 !! or from arbitrary sea-ice conditions 65 !!------------------------------------------------------------------- 66 INTEGER :: ji, jj, jk, jl ! dummy loop indices 67 REAL(wp) :: zeps6, zeps, ztmelts, epsi06 ! local scalars 68 REAL(wp) :: zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs 69 REAL(wp), POINTER, DIMENSION(:) :: zgfactorn, zhin 70 REAL(wp), POINTER, DIMENSION(:) :: zgfactors, zhis 71 REAL(wp), POINTER, DIMENSION(:,:) :: zidto ! ice indicator 72 !-------------------------------------------------------------------- 73 74 CALL wrk_alloc( jpm, zgfactorn, zgfactors, zhin, zhis ) 67 !! ** Method : 68 !! This routine will put some ice where ocean 69 !! is at the freezing point, then fill in ice 70 !! state variables using prescribed initial 71 !! values in the namelist 72 !! 73 !! ** Steps : 74 !! 1) Read namelist 75 !! 2) Basal temperature; ice and hemisphere masks 76 !! 3) Fill in the ice thickness distribution using gaussian 77 !! 4) Fill in space-dependent arrays for state variables 78 !! 5) Diagnostic arrays 79 !! 6) Lateral boundary conditions 80 !! 81 !! History : 82 !! 2.0 ! 01-04 (C. Ethe, G. Madec) Original code 83 !! 3.0 ! 2007 (M. Vancoppenolle) Rewrite for ice cats 84 !! 4.0 ! 09-11 (M. Vancoppenolle) Enhanced version for ice cats 85 !!-------------------------------------------------------------------- 86 87 !! * Local variables 88 INTEGER :: ji, jj, jk, jl ! dummy loop indices 89 REAL(wp) :: epsi06, epsi20, ztmelts 90 INTEGER :: i_hemis, i_fill, jl0 91 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 92 REAL(wp), POINTER, DIMENSION(:) :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini 93 REAL(wp), POINTER, DIMENSION(:,:) :: zht_i_ini, za_i_ini, zv_i_ini 94 REAL(wp), POINTER, DIMENSION(:,:) :: zidto ! ice indicator 95 INTEGER, POINTER, DIMENSION(:,:) :: zhemis ! hemispheric index 96 !-------------------------------------------------------------------- 97 75 98 CALL wrk_alloc( jpi, jpj, zidto ) 76 77 !-------------------------------------------------------------------- 78 ! 1) Preliminary things 79 !-------------------------------------------------------------------- 80 epsi06 = 1.e-6_wp 99 CALL wrk_alloc( jpi, jpj, zhemis ) 100 CALL wrk_alloc( jpl, 2, zht_i_ini, za_i_ini, zv_i_ini ) 101 CALL wrk_alloc( 2, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 102 103 epsi06 = 1.0e-6 104 epsi20 = 1.0e-20 105 IF(lwp) WRITE(numout,*) 106 IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' 107 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 108 109 !-------------------------------------------------------------------- 110 ! 1) Read namelist 111 !-------------------------------------------------------------------- 81 112 82 113 CALL lim_istate_init ! reading the initials parameters of the ice … … 87 118 88 119 !-------------------------------------------------------------------- 89 ! 2) Ice initialization (hi,hs,frld,t_su,sm_i,t_i,t_s) | 90 !-------------------------------------------------------------------- 91 92 IF(lwp) WRITE(numout,*) 93 IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' 94 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 95 120 ! 2) Basal temperature, ice mask and hemispheric index 121 !-------------------------------------------------------------------- 122 123 ! Basal temperature is set to the freezing point of seawater in Celsius 96 124 t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius] 97 125 98 126 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 99 127 DO ji = 1, jpi 100 IF( tsn(ji,jj,1,jp_tem) - t_bo(ji,jj) >= ttest ) THEN ; zidto(ji,jj) = 0. e0! no ice101 ELSE ; zidto(ji,jj) = 1. e0! ice128 IF( tsn(ji,jj,1,jp_tem) - t_bo(ji,jj) >= ttest ) THEN ; zidto(ji,jj) = 0._wp ! no ice 129 ELSE ; zidto(ji,jj) = 1._wp ! ice 102 130 ENDIF 103 131 END DO 104 132 END DO 105 133 106 t_bo(:,:) = t_bo(:,:) + rt0 ! t_bo converted from Celsius to Kelvin (rt0 over land) 107 108 ! constants for heat contents 109 zeps = 1.e-20_wp 110 zeps6 = 1.e-06_wp 111 112 ! zgfactor for initial ice distribution 113 zgfactorn(:) = 0._wp 114 zgfactors(:) = 0._wp 115 116 ! first ice type 117 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 118 zhin (1) = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 119 zgfactorn(1) = zgfactorn(1) + exp(-(zhin(1)-hginn_u)*(zhin(1)-hginn_u) * 0.5_wp ) 120 zhis (1) = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 121 zgfactors(1) = zgfactors(1) + exp(-(zhis(1)-hgins_u)*(zhis(1)-hgins_u) * 0.5_wp ) 122 END DO ! jl 123 zgfactorn(1) = aginn_u / zgfactorn(1) 124 zgfactors(1) = agins_u / zgfactors(1) 125 126 ! ------------- 127 ! new distribution, polynom of second order, conserving area and volume 128 zh1 = 0._wp 129 zh2 = 0._wp 130 zh3 = 0._wp 131 DO jl = 1, jpl 132 zh = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 133 zh1 = zh1 + zh 134 zh2 = zh2 + zh * zh 135 zh3 = zh3 + zh * zh * zh 136 END DO 137 IF(lwp) WRITE(numout,*) ' zh1 : ', zh1 138 IF(lwp) WRITE(numout,*) ' zh2 : ', zh2 139 IF(lwp) WRITE(numout,*) ' zh3 : ', zh3 140 141 zvol = aginn_u * hginn_u 142 zare = aginn_u 143 IF( jpl >= 2 ) THEN 144 zbn = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3) 145 zan = ( zare - zbn*zh1 ) / zh2 146 ENDIF 147 148 IF(lwp) WRITE(numout,*) ' zvol: ', zvol 149 IF(lwp) WRITE(numout,*) ' zare: ', zare 150 IF(lwp) WRITE(numout,*) ' zbn : ', zbn 151 IF(lwp) WRITE(numout,*) ' zan : ', zan 152 153 zvol = agins_u * hgins_u 154 zare = agins_u 155 IF( jpl >= 2 ) THEN 156 zbs = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3) 157 zas = ( zare - zbs*zh1 ) / zh2 158 ENDIF 159 160 IF(lwp) WRITE(numout,*) ' zvol: ', zvol 161 IF(lwp) WRITE(numout,*) ' zare: ', zare 162 IF(lwp) WRITE(numout,*) ' zbn : ', zbn 163 IF(lwp) WRITE(numout,*) ' zan : ', zan 164 165 !end of new lines 166 ! ------------- 167 !!! 168 ! retour a LIMA_MEC 169 ! ! second ice type 170 ! zdummy = hi_max(ice_cat_bounds(2,1)-1) 171 ! hi_max(ice_cat_bounds(2,1)-1) = 0.0 172 173 ! ! here to change !!!! 174 ! jm = 2 175 ! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 176 ! zhin (2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 177 ! zhin (2) = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm ) + & 178 ! hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm) ) / 2.0 179 ! zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0) 180 ! zhis (2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 181 ! zhis (2) = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm ) + & 182 ! hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm) ) / 2.0 183 ! zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0) 184 ! END DO ! jl 185 ! zgfactorn(2) = aginn_d / zgfactorn(2) 186 ! zgfactors(2) = agins_d / zgfactors(2) 187 ! hi_max(ice_cat_bounds(2,1)-1) = zdummy 188 ! END retour a LIMA_MEC 189 !!! 190 191 !!gm optimisation : loop over the ice categories inside the ji, jj loop !!! 192 134 t_bo(:,:) = t_bo(:,:) + rt0 ! conversion to Kelvin 135 136 ! Hemispheric index 137 ! MV 2011 new initialization 193 138 DO jj = 1, jpj 194 139 DO ji = 1, jpi 195 196 !--- Northern hemisphere197 !----------------------------------------------------------------198 140 IF( fcor(ji,jj) >= 0._wp ) THEN 199 200 !----------------------- 201 ! Ice area / thickness 202 !----------------------- 203 204 IF ( jpl .EQ. 1) THEN ! one category 205 206 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 207 a_i(ji,jj,jl) = zidto(ji,jj) * aginn_u 208 ht_i(ji,jj,jl) = zidto(ji,jj) * hginn_u 209 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 210 END DO 211 212 ELSE ! several categories 213 214 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 215 zhin(1) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 216 a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* & 217 (zhin(1)-hginn_u)/2.0) , epsi06) 218 ! new line 219 a_i(ji,jj,jl) = zidto(ji,jj) * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) ) 220 ht_i(ji,jj,jl) = zidto(ji,jj) * zhin(1) 221 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 222 END DO 223 224 ENDIF 225 226 227 !!! 228 ! retour a LIMA_MEC 229 ! !ridged ice 230 ! zdummy = hi_max(ice_cat_bounds(2,1)-1) 231 ! hi_max(ice_cat_bounds(2,1)-1) = 0.0 232 ! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories 233 ! zhin(2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 234 ! a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* & 235 ! (zhin(2)-hginn_d)/2.0) , epsi06) 236 ! ht_i(ji,jj,jl) = zidto(ji,jj) * zhin(2) 237 ! v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 238 ! END DO 239 ! hi_max(ice_cat_bounds(2,1)-1) = zdummy 240 241 ! !rafted ice 242 ! jl = 6 243 ! a_i(ji,jj,jl) = 0.0 244 ! ht_i(ji,jj,jl) = 0.0 245 ! v_i(ji,jj,jl) = 0.0 246 ! END retour a LIMA_MEC 247 !!! 248 249 DO jl = 1, jpl 250 251 !------------- 252 ! Snow depth 253 !------------- 254 ht_s(ji,jj,jl) = zidto(ji,jj) * hninn 255 v_s(ji,jj,jl) = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 256 257 !--------------- 258 ! Ice salinity 259 !--------------- 260 sm_i(ji,jj,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1 261 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 262 263 !---------- 264 ! Ice age 265 !---------- 266 o_i(ji,jj,jl) = zidto(ji,jj) * 1.0 + ( 1.0 - zidto(ji,jj) ) 267 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) 268 269 !------------------------------ 270 ! Sea ice surface temperature 271 !------------------------------ 272 273 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 274 275 !------------------------------------ 276 ! Snow temperature and heat content 277 !------------------------------------ 278 279 DO jk = 1, nlay_s 280 t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 281 ! Snow energy of melting 282 e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 283 ! Change dimensions 284 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 285 ! Multiply by volume, so that heat content in 10^9 Joules 286 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 287 v_s(ji,jj,jl) / nlay_s 288 END DO !jk 289 290 !----------------------------------------------- 291 ! Ice salinities, temperature and heat content 292 !----------------------------------------------- 293 294 DO jk = 1, nlay_i 295 t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 296 s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1 297 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 298 299 ! heat content per unit volume 300 e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 301 ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 302 + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 303 - rcp * ( ztmelts - rtt ) & 304 ) 305 306 ! Correct dimensions to avoid big values 307 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 308 309 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 310 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 311 area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 312 nlay_i 313 END DO ! jk 314 315 END DO ! jl 316 317 ELSE ! on fcor 318 319 !--- Southern hemisphere 320 !---------------------------------------------------------------- 321 322 !----------------------- 323 ! Ice area / thickness 324 !----------------------- 325 326 IF ( jpl .EQ. 1) THEN ! one category 327 328 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 329 a_i(ji,jj,jl) = zidto(ji,jj) * agins_u 330 ht_i(ji,jj,jl) = zidto(ji,jj) * hgins_u 331 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 332 END DO 333 334 ELSE ! several categories 335 336 !level ice 337 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories 338 339 zhis(1) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 340 a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * & 341 (zhis(1)-hgins_u)/2.0) , epsi06 ) 342 ! new line square distribution volume conserving 343 a_i(ji,jj,jl) = zidto(ji,jj) * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) ) 344 ht_i(ji,jj,jl) = zidto(ji,jj) * zhis(1) 345 v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 346 347 END DO ! jl 348 349 ENDIF 350 351 !!! 352 ! retour a LIMA_MEC 353 ! !ridged ice 354 ! zdummy = hi_max(ice_cat_bounds(2,1)-1) 355 ! hi_max(ice_cat_bounds(2,1)-1) = 0.0 356 ! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories 357 ! zhis(2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 358 ! a_i(ji,jj,jl) = zidto(ji,jj)*MAX( zgfactors(2) & 359 ! & * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 ) 360 ! ht_i(ji,jj,jl) = zidto(ji,jj) * zhis(2) 361 ! v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 362 ! END DO 363 ! hi_max(ice_cat_bounds(2,1)-1) = zdummy 364 365 ! !rafted ice 366 ! jl = 6 367 ! a_i(ji,jj,jl) = 0.0 368 ! ht_i(ji,jj,jl) = 0.0 369 ! v_i(ji,jj,jl) = 0.0 370 ! END retour a LIMA_MEC 371 !!! 372 373 DO jl = 1, jpl !over thickness categories 374 375 !--------------- 376 ! Snow depth 377 !--------------- 378 379 ht_s(ji,jj,jl) = zidto(ji,jj) * hnins 380 v_s(ji,jj,jl) = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 381 382 !--------------- 383 ! Ice salinity 384 !--------------- 385 386 sm_i(ji,jj,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1 387 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 388 389 !---------- 390 ! Ice age 391 !---------- 392 393 o_i(ji,jj,jl) = zidto(ji,jj) * 1.0 + ( 1.0 - zidto(ji,jj) ) 394 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) 395 396 !------------------------------ 397 ! Sea ice surface temperature 398 !------------------------------ 399 400 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 401 402 !---------------------------------- 403 ! Snow temperature / heat content 404 !---------------------------------- 405 406 DO jk = 1, nlay_s 407 t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 408 ! Snow energy of melting 409 e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 410 ! Change dimensions 411 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 412 ! Multiply by volume, so that heat content in 10^9 Joules 413 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 414 v_s(ji,jj,jl) / nlay_s 415 END DO 416 417 !--------------------------------------------- 418 ! Ice temperature, salinity and heat content 419 !--------------------------------------------- 420 421 DO jk = 1, nlay_i 422 t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 423 s_i(ji,jj,jk,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1 424 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 425 426 ! heat content per unit volume 427 e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 428 ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 429 + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 430 - rcp * ( ztmelts - rtt ) & 431 ) 432 433 ! Correct dimensions to avoid big values 434 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 435 436 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 437 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 438 area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 439 nlay_i 440 END DO !jk 441 442 END DO ! jl 443 444 ENDIF ! on fcor 445 141 zhemis(ji,jj) = 1 ! Northern hemisphere 142 ELSE 143 zhemis(ji,jj) = 2 ! Southern hemisphere 144 ENDIF 446 145 END DO 447 146 END DO 448 449 !-------------------------------------------------------------------- 450 ! 3) Global ice variables for output diagnostics | 451 !-------------------------------------------------------------------- 452 453 fsbbq (:,:) = 0.e0 454 u_ice (:,:) = 0.e0 455 v_ice (:,:) = 0.e0 456 stress1_i(:,:) = 0.0 457 stress2_i(:,:) = 0.0 458 stress12_i(:,:) = 0.0 459 460 !-------------------------------------------------------------------- 461 ! 4) Moments for advection 462 !-------------------------------------------------------------------- 463 464 sxopw (:,:) = 0.e0 465 syopw (:,:) = 0.e0 466 sxxopw(:,:) = 0.e0 467 syyopw(:,:) = 0.e0 468 sxyopw(:,:) = 0.e0 469 470 sxice (:,:,:) = 0.e0 ; sxsn (:,:,:) = 0.e0 ; sxa (:,:,:) = 0.e0 471 syice (:,:,:) = 0.e0 ; sysn (:,:,:) = 0.e0 ; sya (:,:,:) = 0.e0 472 sxxice(:,:,:) = 0.e0 ; sxxsn(:,:,:) = 0.e0 ; sxxa (:,:,:) = 0.e0 473 syyice(:,:,:) = 0.e0 ; syysn(:,:,:) = 0.e0 ; syya (:,:,:) = 0.e0 474 sxyice(:,:,:) = 0.e0 ; sxysn(:,:,:) = 0.e0 ; sxya (:,:,:) = 0.e0 475 476 sxc0 (:,:,:) = 0.e0 ; sxe (:,:,:,:)= 0.e0 477 syc0 (:,:,:) = 0.e0 ; sye (:,:,:,:)= 0.e0 478 sxxc0 (:,:,:) = 0.e0 ; sxxe (:,:,:,:)= 0.e0 479 syyc0 (:,:,:) = 0.e0 ; syye (:,:,:,:)= 0.e0 480 sxyc0 (:,:,:) = 0.e0 ; sxye (:,:,:,:)= 0.e0 481 482 sxsal (:,:,:) = 0.e0 483 sysal (:,:,:) = 0.e0 484 sxxsal (:,:,:) = 0.e0 485 syysal (:,:,:) = 0.e0 486 sxysal (:,:,:) = 0.e0 487 488 !-------------------------------------------------------------------- 489 ! 5) Lateral boundary conditions | 147 ! END MV 2011 new initialization 148 149 !-------------------------------------------------------------------- 150 ! 3) Initialization of sea ice state variables 151 !-------------------------------------------------------------------- 152 153 !----------------------------- 154 ! 3.1) Hemisphere-dependent arrays 155 !----------------------------- 156 ! assign initial thickness, concentration, snow depth and salinity to 157 ! an hemisphere-dependent array 158 zhm_i_ini(1) = hginn ; zhm_i_ini(2) = hgins ! ice thickness 159 zat_i_ini(1) = aginn ; zat_i_ini(2) = agins ! ice concentration 160 zvt_i_ini(:) = zhm_i_ini(:) * zat_i_ini(:) ! ice volume 161 zhm_s_ini(1) = hninn ; zhm_s_ini(2) = hnins ! snow depth 162 zsm_i_ini(1) = sinn ; zsm_i_ini(2) = sins ! bulk ice salinity 163 164 !--------------------------------------------------------------------- 165 ! 3.2) Distribute ice concentration and thickness into the categories 166 !--------------------------------------------------------------------- 167 ! a gaussian distribution for ice concentration is used 168 ! then we check whether the distribution fullfills 169 ! volume and area conservation, positivity and ice categories bounds 170 DO i_hemis = 1, 2 171 172 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 173 174 ! note for the great nemo engineers: 175 ! only very few of the WRITE statements are necessary for the reference version 176 ! they were one day useful, but now i personally doubt of their 177 ! potential for bringing anything useful 178 179 DO i_fill = jpl, 1, -1 180 IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 181 !---------------------------- 182 ! fill the i_fill categories 183 !---------------------------- 184 ! *** 1 category to fill 185 IF ( i_fill .EQ. 1 ) THEN 186 zht_i_ini(1,i_hemis) = zhm_i_ini(i_hemis) 187 za_i_ini(1,i_hemis) = zat_i_ini(i_hemis) 188 zht_i_ini(2:jpl,i_hemis) = 0._wp 189 za_i_ini(2:jpl,i_hemis) = 0._wp 190 ELSE 191 192 ! *** >1 categores to fill 193 !--- Ice thicknesses in the i_fill - 1 first categories 194 DO jl = 1, i_fill - 1 195 zht_i_ini(jl,i_hemis) = 0.5 * ( hi_max(jl) + hi_max(jl-1) ) 196 END DO 197 198 !--- jl0: most likely index where cc will be maximum 199 DO jl = 1, jpl 200 IF ( ( zhm_i_ini(i_hemis) .GT. hi_max(jl-1) ) .AND. & 201 ( zhm_i_ini(i_hemis) .LE. hi_max(jl) ) ) THEN 202 jl0 = jl 203 ENDIF 204 END DO 205 jl0 = MIN(jl0, i_fill) 206 207 !--- Concentrations 208 za_i_ini(jl0,i_hemis) = zat_i_ini(i_hemis) / SQRT(REAL(jpl)) 209 DO jl = 1, i_fill - 1 210 IF ( jl .NE. jl0 ) THEN 211 zsigma = 0.5 * zhm_i_ini(i_hemis) 212 zarg = ( zht_i_ini(jl,i_hemis) - zhm_i_ini(i_hemis) ) / zsigma 213 za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2) 214 ENDIF 215 END DO 216 217 zA = 0. ! sum of the areas in the jpl categories 218 DO jl = 1, i_fill - 1 219 zA = zA + za_i_ini(jl,i_hemis) 220 END DO 221 za_i_ini(i_fill,i_hemis) = zat_i_ini(i_hemis) - zA ! ice conc in the last category 222 IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 223 224 !--- Ice thickness in the last category 225 zV = 0. ! sum of the volumes of the N-1 categories 226 DO jl = 1, i_fill - 1 227 zV = zV + za_i_ini(jl,i_hemis)*zht_i_ini(jl,i_hemis) 228 END DO 229 zht_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 230 IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 231 232 !--- volumes 233 zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zht_i_ini(:,i_hemis) 234 IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 235 236 ENDIF ! i_fill 237 238 !--------------------- 239 ! Compatibility tests 240 !--------------------- 241 ! Test 1: area conservation 242 zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons ) 243 IF ( zconv .LT. 1.0e-6 ) THEN 244 ztest_1 = 1 245 ELSE 246 ! this write is useful 247 IF(lwp) WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 248 ztest_1 = 0 249 ENDIF 250 251 ! Test 2: volume conservation 252 zV_cons = SUM(zv_i_ini(:,i_hemis)) 253 zconv = ABS(zvt_i_ini(i_hemis) - zV_cons) 254 255 IF ( zconv .LT. 1.0e-6 ) THEN 256 ztest_2 = 1 257 ELSE 258 ! this write is useful 259 IF(lwp) WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 260 ' zvt_i_ini = ', zvt_i_ini(i_hemis) 261 ztest_2 = 0 262 ENDIF 263 264 ! Test 3: thickness of the last category is in-bounds ? 265 IF ( zht_i_ini(i_fill, i_hemis) .GT. hi_max(i_fill-1) ) THEN 266 ztest_3 = 1 267 ELSE 268 ! this write is useful 269 IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zht_i_ini(i_fill,i_hemis) = ', & 270 zht_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 271 ztest_3 = 0 272 ENDIF 273 274 ! Test 4: positivity of ice concentrations 275 ztest_4 = 1 276 DO jl = 1, jpl 277 IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 278 ! this write is useful 279 IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis) 280 ztest_4 = 0 281 ENDIF 282 END DO 283 284 ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 285 286 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 287 288 END DO ! i_fill 289 290 IF(lwp) THEN 291 WRITE(numout,*), ' ztests : ', ztests 292 IF ( ztests .NE. 4 ) THEN 293 WRITE(numout,*) 294 WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 295 WRITE(numout,*), ' !!!! RED ALERT !!! ' 296 WRITE(numout,*), ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!' 297 WRITE(numout,*), ' !!!! Something is wrong in the LIM3 initialization procedure ' 298 WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 299 WRITE(numout,*) 300 WRITE(numout,*), ' *** ztests is not equal to 4 ' 301 WRITE(numout,*), ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 302 WRITE(numout,*), ' zat_i_ini : ', zat_i_ini(i_hemis) 303 WRITE(numout,*), ' zhm_i_ini : ', zhm_i_ini(i_hemis) 304 ENDIF ! ztests .NE. 4 305 ENDIF 306 307 END DO ! i_hemis 308 309 !--------------------------------------------------------------------- 310 ! 3.3) Space-dependent arrays for ice state variables 311 !--------------------------------------------------------------------- 312 313 ! Ice concentration, thickness and volume, snow depth, ice 314 ! salinity, ice age, surface temperature 315 DO jl = 1, jpl ! loop over categories 316 DO jj = 1, jpj 317 DO ji = 1, jpi 318 a_i(ji,jj,jl) = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj)) ! concentration 319 ht_i(ji,jj,jl) = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj)) ! ice thickness 320 ht_s(ji,jj,jl) = zidto(ji,jj) * zhm_s_ini(zhemis(ji,jj)) ! snow depth 321 sm_i(ji,jj,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) ! salinity 322 o_i(ji,jj,jl) = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age 323 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * t_bo(ji,jj) ! surf temp 324 325 ! ice volume, snow volume, salt content, age content 326 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume 327 v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume 328 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 329 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content 330 END DO ! ji 331 END DO ! jj 332 END DO ! jl 333 334 ! Snow temperature and heat content 335 DO jk = 1, nlay_s 336 DO jl = 1, jpl ! loop over categories 337 DO jj = 1, jpj 338 DO ji = 1, jpi 339 t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * rtt 340 ! Snow energy of melting 341 e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 342 ! Change dimensions 343 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 344 ! Multiply by volume, so that heat content in 10^9 Joules 345 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 346 END DO ! ji 347 END DO ! jj 348 END DO ! jl 349 END DO ! jk 350 351 ! Ice salinity, temperature and heat content 352 DO jk = 1, nlay_i 353 DO jl = 1, jpl ! loop over categories 354 DO jj = 1, jpj 355 DO ji = 1, jpi 356 t_i(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1._wp - zidto(ji,jj) ) * rtt 357 s_i(ji,jj,jk,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min 358 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 359 360 ! heat content per unit volume 361 e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 362 + lfus * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 363 - rcp * ( ztmelts - rtt ) ) 364 365 ! Correct dimensions to avoid big values 366 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 367 368 ! Mutliply by ice volume, and divide by number of layers 369 ! to get heat content in 10^9 J 370 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i 371 END DO ! ji 372 END DO ! jj 373 END DO ! jl 374 END DO ! jk 375 376 !-------------------------------------------------------------------- 377 ! 4) Global ice variables for output diagnostics | 378 !-------------------------------------------------------------------- 379 fsbbq (:,:) = 0._wp 380 u_ice (:,:) = 0._wp 381 v_ice (:,:) = 0._wp 382 stress1_i(:,:) = 0._wp 383 stress2_i(:,:) = 0._wp 384 stress12_i(:,:) = 0._wp 385 386 # if defined key_coupled 387 albege(:,:) = 0.8 * tms(:,:) 388 # endif 389 390 !-------------------------------------------------------------------- 391 ! 5) Moments for advection 392 !-------------------------------------------------------------------- 393 394 sxopw (:,:) = 0._wp 395 syopw (:,:) = 0._wp 396 sxxopw(:,:) = 0._wp 397 syyopw(:,:) = 0._wp 398 sxyopw(:,:) = 0._wp 399 400 sxice (:,:,:) = 0._wp ; sxsn (:,:,:) = 0._wp ; sxa (:,:,:) = 0._wp 401 syice (:,:,:) = 0._wp ; sysn (:,:,:) = 0._wp ; sya (:,:,:) = 0._wp 402 sxxice(:,:,:) = 0._wp ; sxxsn(:,:,:) = 0._wp ; sxxa (:,:,:) = 0._wp 403 syyice(:,:,:) = 0._wp ; syysn(:,:,:) = 0._wp ; syya (:,:,:) = 0._wp 404 sxyice(:,:,:) = 0._wp ; sxysn(:,:,:) = 0._wp ; sxya (:,:,:) = 0._wp 405 406 sxc0 (:,:,:) = 0._wp ; sxe (:,:,:,:)= 0._wp 407 syc0 (:,:,:) = 0._wp ; sye (:,:,:,:)= 0._wp 408 sxxc0 (:,:,:) = 0._wp ; sxxe (:,:,:,:)= 0._wp 409 syyc0 (:,:,:) = 0._wp ; syye (:,:,:,:)= 0._wp 410 sxyc0 (:,:,:) = 0._wp ; sxye (:,:,:,:)= 0._wp 411 412 sxsal (:,:,:) = 0._wp 413 sysal (:,:,:) = 0._wp 414 sxxsal (:,:,:) = 0._wp 415 syysal (:,:,:) = 0._wp 416 sxysal (:,:,:) = 0._wp 417 418 !-------------------------------------------------------------------- 419 ! 6) Lateral boundary conditions | 490 420 !-------------------------------------------------------------------- 491 421 492 422 DO jl = 1, jpl 423 493 424 CALL lbc_lnk( a_i(:,:,jl) , 'T', 1. ) 494 425 CALL lbc_lnk( v_i(:,:,jl) , 'T', 1. ) … … 496 427 CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. ) 497 428 CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. ) 498 ! 429 499 430 CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. ) 500 431 CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. ) … … 513 444 a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl) 514 445 END DO 446 447 at_i (:,:) = 0.0_wp 448 DO jl = 1, jpl 449 at_i (:,:) = at_i (:,:) + a_i (:,:,jl) 450 END DO 515 451 516 452 CALL lbc_lnk( at_i , 'T', 1. ) … … 519 455 CALL lbc_lnk( fsbbq , 'T', 1. ) 520 456 ! 521 CALL wrk_dealloc( jpm, zgfactorn, zgfactors, zhin, zhis ) 457 !-------------------------------------------------------------------- 458 ! 6) ???? | 459 !-------------------------------------------------------------------- 460 tn_ice (:,:,:) = t_su (:,:,:) 461 522 462 CALL wrk_dealloc( jpi, jpj, zidto ) 523 ! 463 CALL wrk_dealloc( jpi, jpj, zhemis ) 464 CALL wrk_dealloc( jpl, 2, zht_i_ini, za_i_ini, zv_i_ini ) 465 CALL wrk_dealloc( 2, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 466 524 467 END SUBROUTINE lim_istate 525 526 468 527 469 SUBROUTINE lim_istate_init … … 531 473 !! ** Purpose : Definition of initial state of the ice 532 474 !! 533 !! ** Method : Read the namiceini namelist and check the parameter 534 !! values called at the first timestep (nit000) 535 !! 536 !! ** input : namelist namiceini 475 !! ** Method : Read the namiceini namelist and check the parameter 476 !! values called at the first timestep (nit000) 477 !! 478 !! ** input : 479 !! Namelist namiceini 480 !! 481 !! history : 482 !! 8.5 ! 03-08 (C. Ethe) original code 483 !! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization 537 484 !!----------------------------------------------------------------------------- 485 <<<<<<< .courant 538 486 INTEGER :: ios ! Local integer output status for namelist read 539 487 NAMELIST/namiceini/ ttest, hninn, hginn_u, aginn_u, hginn_d, aginn_d, hnins, & 540 488 & hgins_u, agins_u, hgins_d, agins_d, sinn, sins 489 ======= 490 NAMELIST/namiceini/ ttest, hninn, hnins, hginn, hgins, aginn, agins, sinn, sins 491 >>>>>>> .fusion-droit.r4160 541 492 !!----------------------------------------------------------------------------- 493 <<<<<<< .courant 542 494 ! 543 495 REWIND( numnam_ice_ref ) ! Namelist namiceini in reference namelist : Ice initial state … … 549 501 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp ) 550 502 WRITE ( numoni, namiceini ) 551 ! 552 IF(lwp) THEN ! control print 503 ======= 504 505 ! Define the initial parameters 506 ! ------------------------- 507 508 ! Read Namelist namiceini 509 REWIND ( numnam_ice ) 510 READ ( numnam_ice , namiceini ) 511 >>>>>>> .fusion-droit.r4160 512 IF(lwp) THEN 553 513 WRITE(numout,*) 554 514 WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' … … 556 516 WRITE(numout,*) ' threshold water temp. for initial sea-ice ttest = ', ttest 557 517 WRITE(numout,*) ' initial snow thickness in the north hninn = ', hninn 558 WRITE(numout,*) ' initial undef ice thickness in the north hginn_u = ', hginn_u559 WRITE(numout,*) ' initial undef ice concentr. in the north aginn_u = ', aginn_u560 WRITE(numout,*) ' initial def ice thickness in the north hginn_d = ', hginn_d561 WRITE(numout,*) ' initial def ice concentr. in the north aginn_d = ', aginn_d562 518 WRITE(numout,*) ' initial snow thickness in the south hnins = ', hnins 563 WRITE(numout,*) ' initial undef ice thickness in the north hgins_u = ', hgins_u564 WRITE(numout,*) ' initial undef ice concentr. in the north agins_u = ', agins_u565 WRITE(numout,*) ' initial def ice thickness in the north hgins_d = ', hgins_d566 WRITE(numout,*) ' initial def ice concentr. in the north agins_d = ', agins_d567 WRITE(numout,*) ' initial ice salinity in the northsinn = ', sinn568 WRITE(numout,*) ' initial ice salinity in the southsins = ', sins519 WRITE(numout,*) ' initial ice thickness in the north hginn = ', hginn 520 WRITE(numout,*) ' initial ice thickness in the south hgins = ', hgins 521 WRITE(numout,*) ' initial ice concentr. in the north aginn = ', aginn 522 WRITE(numout,*) ' initial ice concentr. in the north agins = ', agins 523 WRITE(numout,*) ' initial ice salinity in the north sinn = ', sinn 524 WRITE(numout,*) ' initial ice salinity in the south sins = ', sins 569 525 ENDIF 570 ! 526 571 527 END SUBROUTINE lim_istate_init 572 528 -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r4147 r4161 5 5 !!====================================================================== 6 6 !! History : LIM ! 2006-02 (M. Vancoppenolle) Original code 7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & fsalt_rpo7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_mec 8 8 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 9 9 !!---------------------------------------------------------------------- … … 12 12 !! 'key_lim3' LIM-3 sea-ice model 13 13 !!---------------------------------------------------------------------- 14 USE par_oce ! ocean parameters 15 USE dom_oce ! ocean domain 16 USE phycst ! physical constants (ocean directory) 17 USE sbc_oce ! surface boundary condition: ocean fields 18 USE thd_ice ! LIM thermodynamics 19 USE ice ! LIM variables 20 USE par_ice ! LIM parameters 21 USE dom_ice ! LIM domain 22 USE limthd_lac ! LIM 23 USE limvar ! LIM 24 USE limcons ! LIM 25 USE in_out_manager ! I/O manager 26 USE lbclnk ! lateral boundary condition - MPP exchanges 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! work arrays 29 USE prtctl ! Print control 30 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 14 USE par_oce ! ocean parameters 15 USE dom_oce ! ocean domain 16 USE phycst ! physical constants (ocean directory) 17 USE sbc_oce ! surface boundary condition: ocean fields 18 USE thd_ice ! LIM thermodynamics 19 USE ice ! LIM variables 20 USE par_ice ! LIM parameters 21 USE dom_ice ! LIM domain 22 USE limthd_lac ! LIM 23 USE limvar ! LIM 24 USE limcons ! LIM 25 USE in_out_manager ! I/O manager 26 USE lbclnk ! lateral boundary condition - MPP exchanges 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! work arrays 29 USE prtctl ! Print control 30 ! Check budget (Rousset) 31 USE iom ! I/O manager 32 USE lib_fortran ! glob_sum 33 USE limdiahsb 34 USE timing ! Timing 32 35 33 36 IMPLICIT NONE … … 62 65 REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier 63 66 REAL(wp), PARAMETER :: kraft = 2.0_wp ! rafting multipliyer 67 REAL(wp), PARAMETER :: kamax = 1.0 64 68 65 69 REAL(wp) :: Cp ! … … 141 145 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2) 142 146 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 147 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 148 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 149 ! mass and salt flux (clem) 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... 143 151 !!----------------------------------------------------------------------------- 152 IF( nn_timing == 1 ) CALL timing_start('limitd_me') 144 153 145 154 CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 155 156 CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem 146 157 147 158 IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only) … … 151 162 CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ') 152 163 ENDIF 164 165 IF( ln_limdyn ) THEN ! Start ridging and rafting ! 166 ! ------------------------------- 167 !- check conservation (C Rousset) 168 IF (ln_limdiahsb) THEN 169 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 170 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 171 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 172 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 173 ENDIF 174 !- check conservation (C Rousset) 175 ! ------------------------------- 176 177 ! mass and salt flux init (clem) 178 zviold(:,:,:) = v_i(:,:,:) 179 zvsold(:,:,:) = v_s(:,:,:) 180 zsmvold(:,:,:) = smv_i(:,:,:) 153 181 154 182 !-----------------------------------------------------------------------------! … … 204 232 ! to give asum = 1.0 after ridging. 205 233 206 divu_adv(ji,jj) = ( 1._wp- asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep234 divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep 207 235 208 236 IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) … … 288 316 DO jj = 1, jpj 289 317 DO ji = 1, jpi 290 IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN318 IF (ABS(asum(ji,jj) - kamax ) .LT. epsi11) THEN 291 319 closing_net(ji,jj) = 0._wp 292 320 opning (ji,jj) = 0._wp 293 321 ELSE 294 322 iterate_ridging = 1 295 divu_adv (ji,jj) = ( 1._wp - asum(ji,jj)) * r1_rdtice323 divu_adv (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice 296 324 closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 297 325 opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) ) … … 330 358 DO ji = 1, jpi 331 359 332 IF( ABS( asum(ji,jj) - 1.0 ) > epsi11 )asum_error = .true.360 IF(ABS(asum(ji,jj) - kamax) > epsi11 ) asum_error = .true. 333 361 334 362 dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice … … 349 377 DO jj = 1, jpj 350 378 DO ji = 1, jpi 351 IF( ABS( asum(ji,jj) - 1._wp) > epsi11 ) THEN ! there is a bug379 IF( ABS( asum(ji,jj) - kamax) > epsi11 ) THEN ! there is a bug 352 380 WRITE(numout,*) ' ' 353 381 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) … … 377 405 CALL lim_var_glo2eqv 378 406 CALL lim_itd_me_zapsmall 407 408 !-------------------------------- 409 ! Update mass/salt fluxes (clem) 410 !-------------------------------- 411 DO jl = 1, jpl 412 DO jj = 1, jpj 413 DO ji = 1, jpi 414 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 415 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 416 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 417 sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice 418 END DO 419 END DO 420 END DO 379 421 380 422 !----------------- … … 425 467 ENDIF 426 468 469 ! ------------------------------- 470 !- check conservation (C Rousset) 471 IF (ln_limdiahsb) THEN 472 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 473 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 474 475 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 476 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 477 478 zchk_vmin = glob_min(v_i) 479 zchk_amax = glob_max(SUM(a_i,dim=3)) 480 zchk_amin = glob_min(a_i) 481 482 IF(lwp) THEN 483 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_me) = ',(zchk_v_i * rday) 484 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday) 485 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_me) = ',(zchk_vmin * 1.e-3) 486 IF ( zchk_amax > kamax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_me) = ',zchk_amax 487 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_me) = ',zchk_amin 488 ENDIF 489 ENDIF 490 !- check conservation (C Rousset) 491 ! ------------------------------- 492 427 493 !-------------------------! 428 494 ! Back to initial values … … 448 514 449 515 ! heat content has to be corrected before ice volume 450 DO jl = 1, jpl 451 DO jk = 1, nlay_i 452 DO jj = 1, jpj 453 DO ji = 1, jpi 454 IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 455 ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 456 old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl) 457 d_e_i_trp(ji,jj,jk,jl) = 0._wp 458 ENDIF 459 END DO 460 END DO 461 END DO 462 END DO 463 464 DO jl = 1, jpl 465 DO jj = 1, jpj 466 DO ji = 1, jpi 467 IF( old_v_i (ji,jj,jl) < epsi06 .AND. & 468 d_v_i_trp(ji,jj,jl) > epsi06 ) THEN 469 old_v_i (ji,jj,jl) = d_v_i_trp(ji,jj,jl) 470 d_v_i_trp (ji,jj,jl) = 0._wp 471 old_a_i (ji,jj,jl) = d_a_i_trp(ji,jj,jl) 472 d_a_i_trp (ji,jj,jl) = 0._wp 473 old_v_s (ji,jj,jl) = d_v_s_trp(ji,jj,jl) 474 d_v_s_trp (ji,jj,jl) = 0._wp 475 old_e_s (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 476 d_e_s_trp (ji,jj,1,jl) = 0._wp 477 old_oa_i (ji,jj,jl) = d_oa_i_trp(ji,jj,jl) 478 d_oa_i_trp(ji,jj,jl) = 0._wp 479 IF( num_sal == 2 ) old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 480 d_smv_i_trp(ji,jj,jl) = 0._wp 481 ENDIF 482 END DO 483 END DO 484 END DO 516 !clem@order 517 ! DO jl = 1, jpl 518 ! DO jk = 1, nlay_i 519 ! DO jj = 1, jpj 520 ! DO ji = 1, jpi 521 ! IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 522 ! ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 523 ! old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl) 524 ! d_e_i_trp(ji,jj,jk,jl) = 0._wp 525 ! ENDIF 526 ! END DO 527 ! END DO 528 ! END DO 529 ! END DO 530 ! 531 ! DO jl = 1, jpl 532 ! DO jj = 1, jpj 533 ! DO ji = 1, jpi 534 ! IF( old_v_i (ji,jj,jl) < epsi06 .AND. & 535 ! d_v_i_trp(ji,jj,jl) > epsi06 ) THEN 536 ! old_v_i (ji,jj,jl) = d_v_i_trp(ji,jj,jl) 537 ! d_v_i_trp (ji,jj,jl) = 0._wp 538 ! old_a_i (ji,jj,jl) = d_a_i_trp(ji,jj,jl) 539 ! d_a_i_trp (ji,jj,jl) = 0._wp 540 ! old_v_s (ji,jj,jl) = d_v_s_trp(ji,jj,jl) 541 ! d_v_s_trp (ji,jj,jl) = 0._wp 542 ! old_e_s (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 543 ! d_e_s_trp (ji,jj,1,jl) = 0._wp 544 ! old_oa_i (ji,jj,jl) = d_oa_i_trp(ji,jj,jl) 545 ! d_oa_i_trp(ji,jj,jl) = 0._wp 546 ! IF( num_sal == 2 ) old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 547 ! d_smv_i_trp(ji,jj,jl) = 0._wp 548 ! ENDIF 549 ! END DO 550 ! END DO 551 ! END DO 552 !clem@order 553 ENDIF ! ln_limdyn=.true. 485 554 ! 486 555 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 487 556 ! 557 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem 558 ! 559 IF( nn_timing == 1 ) CALL timing_stop('limitd_me') 488 560 END SUBROUTINE lim_itd_me 489 561 … … 1086 1158 afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 1087 1159 1088 IF (afrac(ji,jj) > 1.0+ epsi11) THEN !riging1160 IF (afrac(ji,jj) > kamax + epsi11) THEN !riging 1089 1161 large_afrac = .true. 1090 ELSEIF (afrac(ji,jj) > 1.0) THEN ! roundoff error1091 afrac(ji,jj) = 1.01162 ELSEIF (afrac(ji,jj) > kamax) THEN ! roundoff error 1163 afrac(ji,jj) = kamax 1092 1164 ENDIF 1093 IF (afrft(ji,jj) > 1.0+ epsi11) THEN !rafting1165 IF (afrft(ji,jj) > kamax + epsi11) THEN !rafting 1094 1166 large_afrft = .true. 1095 ELSEIF (afrft(ji,jj) > 1.0) THEN ! roundoff error1096 afrft(ji,jj) = 1.01167 ELSEIF (afrft(ji,jj) > kamax) THEN ! roundoff error 1168 afrft(ji,jj) = kamax 1097 1169 ENDIF 1098 1170 … … 1137 1209 1138 1210 ! ! excess of salt is flushed into the ocean 1139 sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice1140 1141 rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic / rau0 ! increase in ice volume du to seawater frozen in voids1142 1211 !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 1212 1213 !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic ! gurvan: increase in ice volume du to seawater frozen in voids 1214 1143 1215 !------------------------------------ 1144 1216 ! 3.6 Increment ridging diagnostics … … 1150 1222 dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 1151 1223 dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 1152 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice1224 !clem diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 1153 1225 opening (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 1154 1226 … … 1217 1289 1218 1290 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 1219 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / nlay_i1291 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 1220 1292 1221 1293 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) … … 1240 1312 ji = indxi(ij) 1241 1313 jj = indxj(ij) 1242 IF( afrac(ji,jj) > 1.0+ epsi11 ) THEN1314 IF( afrac(ji,jj) > kamax + epsi11 ) THEN 1243 1315 WRITE(numout,*) '' 1244 1316 WRITE(numout,*) ' ardg > a_i' … … 1252 1324 ji = indxi(ij) 1253 1325 jj = indxj(ij) 1254 IF( afrft(ji,jj) > 1.0+ epsi11 ) THEN1326 IF( afrft(ji,jj) > kamax + epsi11 ) THEN 1255 1327 WRITE(numout,*) '' 1256 1328 WRITE(numout,*) ' arft > a_i' -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r3764 r4161 19 19 !! lim_itd_shiftice : 20 20 !!---------------------------------------------------------------------- 21 USE par_oce ! ocean parameters22 USE dom_oce ! ocean domain23 USE phycst ! physical constants (ocean directory)24 USE ice ! LIM-3 variables25 USE par_ice ! LIM-3 parameters26 USE dom_ice ! LIM-3 domain27 USE thd_ice ! LIM-3 thermodynamic variables28 USE limthd_lac ! LIM-3 lateral accretion29 USE limvar ! LIM-3 variables30 USE limcons ! LIM-3 conservation31 USE prtctl ! Print control32 USE in_out_manager ! I/O manager33 USE lib_mpp ! MPP library34 USE wrk_nemo ! work arrays35 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)36 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)21 USE dom_ice ! LIM-3 domain 22 USE par_oce ! ocean parameters 23 USE dom_oce ! ocean domain 24 USE phycst ! physical constants (ocean directory) 25 USE thd_ice ! LIM-3 thermodynamic variables 26 USE ice ! LIM-3 variables 27 USE par_ice ! LIM-3 parameters 28 USE limthd_lac ! LIM-3 lateral accretion 29 USE limvar ! LIM-3 variables 30 USE limcons ! LIM-3 conservation 31 USE prtctl ! Print control 32 USE in_out_manager ! I/O manager 33 USE lib_mpp ! MPP library 34 USE wrk_nemo ! work arrays 35 USE lib_fortran ! to use key_nosignedzero 36 USE timing ! Timing 37 37 38 38 IMPLICIT NONE … … 45 45 PUBLIC lim_itd_shiftice 46 46 47 REAL(wp) :: epsi20 = 1 e-20_wp ! constant values48 REAL(wp) :: epsi13 = 1 e-13_wp !49 REAL(wp) :: epsi10 = 1 e-10_wp !47 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values 48 REAL(wp) :: epsi13 = 1.e-13_wp ! 49 REAL(wp) :: epsi10 = 1.e-10_wp ! 50 50 51 51 !!---------------------------------------------------------------------- 52 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2010)52 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 53 53 !! $Id$ 54 54 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 67 67 ! 68 68 INTEGER :: jl, ja, jm, jbnd1, jbnd2 ! ice types dummy loop index 69 70 !!------------------------------------------------------------------ 69 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 70 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 71 !!------------------------------------------------------------------ 72 IF( nn_timing == 1 ) CALL timing_start('limitd_th') 73 74 ! ------------------------------- 75 !- check conservation (C Rousset) 76 IF (ln_limdiahsb) THEN 77 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 78 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 79 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 80 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 81 ENDIF 82 !- check conservation (C Rousset) 83 ! ------------------------------- 71 84 72 85 IF( kt == nit000 .AND. lwp ) THEN … … 107 120 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 108 121 d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 109 122 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 110 123 d_smv_i_thd(:,:,:) = 0._wp 111 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 124 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 125 126 ! diag only (clem) 127 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 112 128 113 129 IF(ln_ctl) THEN ! Control print … … 142 158 END DO 143 159 ENDIF 144 160 ! 161 ! ------------------------------- 162 !- check conservation (C Rousset) 163 IF( ln_limdiahsb ) THEN 164 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 165 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 166 167 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 168 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 169 170 zchk_vmin = glob_min(v_i) 171 zchk_amax = glob_max(SUM(a_i,dim=3)) 172 zchk_amin = glob_min(a_i) 173 174 IF(lwp) THEN 175 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_th) = ',(zchk_v_i * rday) 176 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * rday) 177 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_th) = ',(zchk_vmin * 1.e-3) 178 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_th) = ',zchk_amax 179 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_th) = ',zchk_amin 180 ENDIF 181 ENDIF 182 !- check conservation (C Rousset) 183 ! ------------------------------- 184 ! 145 185 !- Recover Old values 146 a_i(:,:,:) = old_a_i (:,:,:)147 v_s(:,:,:) = old_v_s (:,:,:)148 v_i(:,:,:) = old_v_i (:,:,:)149 e_s(:,:,:,:) = old_e_s (:,:,:,:)150 e_i(:,:,:,:) = old_e_i (:,:,:,:)151 ! 186 a_i(:,:,:) = old_a_i (:,:,:) 187 v_s(:,:,:) = old_v_s (:,:,:) 188 v_i(:,:,:) = old_v_i (:,:,:) 189 e_s(:,:,:,:) = old_e_s (:,:,:,:) 190 e_i(:,:,:,:) = old_e_i (:,:,:,:) 191 !?? oa_i(:,:,:) = old_oa_i(:,:,:) 152 192 IF( num_sal == 2 ) smv_i(:,:,:) = old_smv_i(:,:,:) 153 193 ! 194 IF( nn_timing == 1 ) CALL timing_stop('limitd_th') 154 195 END SUBROUTINE lim_itd_th 155 196 ! … … 172 213 ! 173 214 INTEGER :: ji, jj, jl ! dummy loop index 174 INTEGER :: zji, zjj, nd ! local integer 215 INTEGER :: ii, ij ! 2D corresponding indices to ji 216 INTEGER :: nd ! local integer 175 217 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 176 REAL(wp) :: zx2, zwk2, zda0, zetamax , zhimin! - -218 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 177 219 REAL(wp) :: zx3, zareamin, zindb ! - - 178 220 CHARACTER (len = 15) :: fieldid … … 210 252 CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 211 253 212 zhimin = 0.1 !minimum ice thickness tolerated by the model213 254 zareamin = epsi10 !minimum area in thickness categories tolerated by the conceptors of the model 214 255 … … 240 281 DO jj = 1, jpj 241 282 DO ji = 1, jpi 242 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl) )) !0 if no ice and 1 if yes283 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 243 284 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi10) * zindb 244 zindb = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl) )) !0 if no ice and 1 if yes285 zindb = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 245 286 zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),epsi10) * zindb 246 287 IF( a_i(ji,jj,jl) > 1e-6 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) … … 285 326 DO jl = klbnd, kubnd - 1 286 327 DO ji = 1, nbrem 287 zji = nind_i(ji)288 zjj = nind_j(ji)328 ii = nind_i(ji) 329 ij = nind_j(ji) 289 330 ! 290 IF ( ( zht_i_o( zji,zjj,jl) .GT.epsi10 ) .AND. &291 ( zht_i_o( zji,zjj,jl+1).GT.epsi10 ) ) THEN331 IF ( ( zht_i_o(ii,ij,jl) .GT.epsi10 ) .AND. & 332 ( zht_i_o(ii,ij,jl+1).GT.epsi10 ) ) THEN 292 333 !interpolate between adjacent category growth rates 293 zslope = ( zdhice( zji,zjj,jl+1) - zdhice(zji,zjj,jl) ) / &294 ( zht_i_o ( zji,zjj,jl+1) - zht_i_o (zji,zjj,jl) )295 zhbnew( zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) + &296 zslope * ( hi_max(jl) - zht_i_o( zji,zjj,jl) )297 ELSEIF (zht_i_o( zji,zjj,jl).gt.epsi10) THEN298 zhbnew( zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl)299 ELSEIF (zht_i_o( zji,zjj,jl+1).gt.epsi10) THEN300 zhbnew( zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl+1)334 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / & 335 ( zht_i_o (ii,ij,jl+1) - zht_i_o (ii,ij,jl) ) 336 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + & 337 zslope * ( hi_max(jl) - zht_i_o(ii,ij,jl) ) 338 ELSEIF (zht_i_o(ii,ij,jl).gt.epsi10) THEN 339 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 340 ELSEIF (zht_i_o(ii,ij,jl+1).gt.epsi10) THEN 341 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 301 342 ELSE 302 zhbnew( zji,zjj,jl) = hi_max(jl)343 zhbnew(ii,ij,jl) = hi_max(jl) 303 344 ENDIF 304 345 END DO … … 307 348 DO ji = 1, nbrem 308 349 ! jl, ji 309 zji = nind_i(ji)310 zjj = nind_j(ji)350 ii = nind_i(ji) 351 ij = nind_j(ji) 311 352 ! jl, ji 312 IF ( ( a_i( zji,zjj,jl) .GT.epsi10) .AND. &313 ( ht_i( zji,zjj,jl).GE. zhbnew(zji,zjj,jl) ) &353 IF ( ( a_i(ii,ij,jl) .GT.epsi10) .AND. & 354 ( ht_i(ii,ij,jl).GE. zhbnew(ii,ij,jl) ) & 314 355 ) THEN 315 zremap_flag( zji,zjj) = 0316 ELSEIF ( ( a_i( zji,zjj,jl+1) .GT. epsi10 ) .AND. &317 ( ht_i( zji,zjj,jl+1).LE. zhbnew(zji,zjj,jl) ) &356 zremap_flag(ii,ij) = 0 357 ELSEIF ( ( a_i(ii,ij,jl+1) .GT. epsi10 ) .AND. & 358 ( ht_i(ii,ij,jl+1).LE. zhbnew(ii,ij,jl) ) & 318 359 ) THEN 319 zremap_flag( zji,zjj) = 0360 zremap_flag(ii,ij) = 0 320 361 ENDIF 321 362 322 363 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 323 364 ! jl, ji 324 IF (zhbnew( zji,zjj,jl).gt.hi_max(jl+1)) THEN325 zremap_flag( zji,zjj) = 0365 IF (zhbnew(ii,ij,jl).gt.hi_max(jl+1)) THEN 366 zremap_flag(ii,ij) = 0 326 367 ENDIF 327 368 ! jl, ji 328 IF (zhbnew( zji,zjj,jl).lt.hi_max(jl-1)) THEN329 zremap_flag( zji,zjj) = 0369 IF (zhbnew(ii,ij,jl).lt.hi_max(jl-1)) THEN 370 zremap_flag(ii,ij) = 0 330 371 ENDIF 331 372 ! jl, ji … … 379 420 !- 7.2 Area lost due to melting of thin ice (first category, klbnd) 380 421 DO ji = 1, nbrem 381 zji = nind_i(ji)382 zjj = nind_j(ji)422 ii = nind_i(ji) 423 ij = nind_j(ji) 383 424 384 425 !ji 385 IF (a_i( zji,zjj,klbnd) .gt. epsi10) THEN386 zdh0 = zdhice( zji,zjj,klbnd) !decrease of ice thickness in the lower category426 IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN 427 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 387 428 ! ji, a_i > epsi10 388 429 IF (zdh0 .lt. 0.0) THEN !remove area from category 1 … … 391 432 392 433 !Integrate g(1) from 0 to dh0 to estimate area melted 393 zetamax = MIN(zdh0,hR( zji,zjj,klbnd)) - hL(zji,zjj,klbnd)434 zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd) 394 435 IF (zetamax.gt.0.0) THEN 395 436 zx1 = zetamax 396 437 zx2 = 0.5 * zetamax*zetamax 397 zda0 = g1( zji,zjj,klbnd) * zx2 + g0(zji,zjj,klbnd) * zx1 !ice area removed438 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 398 439 ! Constrain new thickness <= ht_i 399 zdamax = a_i( zji,zjj,klbnd) * &400 (1.0 - ht_i( zji,zjj,klbnd)/zht_i_o(zji,zjj,klbnd)) ! zdamax > 0440 zdamax = a_i(ii,ij,klbnd) * & 441 (1.0 - ht_i(ii,ij,klbnd)/zht_i_o(ii,ij,klbnd)) ! zdamax > 0 401 442 !ice area lost due to melting of thin ice 402 443 zda0 = MIN(zda0, zdamax) 403 444 404 445 ! Remove area, conserving volume 405 ht_i( zji,zjj,klbnd) = ht_i(zji,zjj,klbnd) &406 * a_i( zji,zjj,klbnd) / ( a_i(zji,zjj,klbnd) - zda0 )407 a_i( zji,zjj,klbnd) = a_i(zji,zjj,klbnd) - zda0408 v_i( zji,zjj,klbnd) = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd)446 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 447 * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 448 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 449 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem@useless ? 409 450 ENDIF ! zetamax > 0 410 451 ! ji, a_i > epsi10 … … 412 453 ELSE ! if ice accretion 413 454 ! ji, a_i > epsi10; zdh0 > 0 414 IF ( ntyp .EQ. 1 ) zhbnew( zji,zjj,klbnd-1) = MIN(zdh0,hi_max(klbnd))455 IF ( ntyp .EQ. 1 ) zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 415 456 ! zhbnew was 0, and is shifted to the right to account for thin ice 416 457 ! growth in openwater (F0 = f1) 417 IF ( ntyp .NE. 1 ) zhbnew( zji,zjj,0) = 0458 IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0 418 459 ! in other types there is 419 460 ! no open water growth (F0 = 0) … … 446 487 447 488 DO ji = 1, nbrem 448 zji = nind_i(ji)449 zjj = nind_j(ji)450 451 IF (zhbnew( zji,zjj,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1489 ii = nind_i(ji) 490 ij = nind_j(ji) 491 492 IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 452 493 453 494 ! left and right integration limits in eta space 454 zvetamin(ji) = MAX(hi_max(jl), hL( zji,zjj,jl)) - hL(zji,zjj,jl)455 zvetamax(ji) = MIN(zhbnew( zji,zjj,jl), hR(zji,zjj,jl)) - hL(zji,zjj,jl)456 zdonor( zji,zjj,jl) = jl495 zvetamin(ji) = MAX(hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl) 496 zvetamax(ji) = MIN(zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl) 497 zdonor(ii,ij,jl) = jl 457 498 458 499 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl … … 460 501 ! left and right integration limits in eta space 461 502 zvetamin(ji) = 0.0 462 zvetamax(ji) = MIN(hi_max(jl), hR( zji,zjj,jl+1)) - hL(zji,zjj,jl+1)463 zdonor( zji,zjj,jl) = jl + 1503 zvetamax(ji) = MIN(hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1) 504 zdonor(ii,ij,jl) = jl + 1 464 505 465 506 ENDIF ! zhbnew(jl) > hi_max(jl) … … 475 516 zwk2 = zwk2 * zetamax 476 517 zx3 = 1.0/3.0 * (zwk2 - zwk1) 477 nd = zdonor( zji,zjj,jl)478 zdaice( zji,zjj,jl) = g1(zji,zjj,nd)*zx2 + g0(zji,zjj,nd)*zx1479 zdvice( zji,zjj,jl) = g1(zji,zjj,nd)*zx3 + g0(zji,zjj,nd)*zx2 + &480 zdaice( zji,zjj,jl)*hL(zji,zjj,nd)518 nd = zdonor(ii,ij,jl) 519 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 520 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + & 521 zdaice(ii,ij,jl)*hL(ii,ij,nd) 481 522 482 523 END DO ! ji … … 493 534 494 535 DO ji = 1, nbrem 495 zji = nind_i(ji) 496 zjj = nind_j(ji) 497 IF ( ( zhimin .GT. 0.0 ) .AND. & 498 ( ( a_i(zji,zjj,1) .GT. epsi10 ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) & 499 ) THEN 500 a_i(zji,zjj,1) = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin 501 ht_i(zji,zjj,1) = zhimin 502 v_i(zji,zjj,1) = a_i(zji,zjj,1)*ht_i(zji,zjj,1) 536 ii = nind_i(ji) 537 ij = nind_j(ji) 538 IF ( ( a_i(ii,ij,1) > epsi10 ) .AND. ( ht_i(ii,ij,1) < hiclim ) ) THEN 539 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim 540 ht_i(ii,ij,1) = hiclim 541 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem@useless 503 542 ENDIF 504 543 END DO !ji … … 625 664 626 665 INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices 627 INTEGER :: zji, zjj ! indices when changing from 2D-1D is done666 INTEGER :: ii, ij ! indices when changing from 2D-1D is done 628 667 629 668 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn … … 759 798 760 799 DO ji = 1, nbrem 761 zji = nind_i(ji)762 zjj = nind_j(ji)763 764 jl1 = zdonor( zji,zjj,jl)765 zindb = MAX( 0.0 , SIGN( 1.0 , v_i( zji,zjj,jl1) - epsi10 ) )766 zworka( zji,zjj) = zdvice(zji,zjj,jl) / MAX(v_i(zji,zjj,jl1),epsi10) * zindb800 ii = nind_i(ji) 801 ij = nind_j(ji) 802 803 jl1 = zdonor(ii,ij,jl) 804 zindb = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 805 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * zindb 767 806 IF( jl1 == jl) THEN ; jl2 = jl1+1 768 807 ELSE ; jl2 = jl … … 773 812 !-------------- 774 813 775 a_i( zji,zjj,jl1) = a_i(zji,zjj,jl1) - zdaice(zji,zjj,jl)776 a_i( zji,zjj,jl2) = a_i(zji,zjj,jl2) + zdaice(zji,zjj,jl)814 a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 815 a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 777 816 778 817 !-------------- … … 780 819 !-------------- 781 820 782 v_i( zji,zjj,jl1) = v_i(zji,zjj,jl1) - zdvice(zji,zjj,jl)783 v_i( zji,zjj,jl2) = v_i(zji,zjj,jl2) + zdvice(zji,zjj,jl)821 v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 822 v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 784 823 785 824 !-------------- … … 787 826 !-------------- 788 827 789 zdvsnow = v_s( zji,zjj,jl1) * zworka(zji,zjj)790 v_s( zji,zjj,jl1) = v_s(zji,zjj,jl1) - zdvsnow791 v_s( zji,zjj,jl2) = v_s(zji,zjj,jl2) + zdvsnow828 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 829 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 830 v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow 792 831 793 832 !-------------------- … … 795 834 !-------------------- 796 835 797 zdesnow = e_s( zji,zjj,1,jl1) * zworka(zji,zjj)798 e_s( zji,zjj,1,jl1) = e_s(zji,zjj,1,jl1) - zdesnow799 e_s( zji,zjj,1,jl2) = e_s(zji,zjj,1,jl2) + zdesnow836 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 837 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow 838 e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow 800 839 801 840 !-------------- … … 803 842 !-------------- 804 843 805 zdo_aice = oa_i( zji,zjj,jl1) * zdaice(zji,zjj,jl)806 oa_i( zji,zjj,jl1) = oa_i(zji,zjj,jl1) - zdo_aice807 oa_i( zji,zjj,jl2) = oa_i(zji,zjj,jl2) + zdo_aice844 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 845 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice 846 oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice 808 847 809 848 !-------------- … … 811 850 !-------------- 812 851 813 zdsm_vice = smv_i( zji,zjj,jl1) * zworka(zji,zjj)814 smv_i( zji,zjj,jl1) = smv_i(zji,zjj,jl1) - zdsm_vice815 smv_i( zji,zjj,jl2) = smv_i(zji,zjj,jl2) + zdsm_vice852 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 853 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice 854 smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice 816 855 817 856 !--------------------- … … 819 858 !--------------------- 820 859 821 zdaTsf = t_su( zji,zjj,jl1) * zdaice(zji,zjj,jl)822 zaTsfn( zji,zjj,jl1) = zaTsfn(zji,zjj,jl1) - zdaTsf823 zaTsfn( zji,zjj,jl2) = zaTsfn(zji,zjj,jl2) + zdaTsf860 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 861 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf 862 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 824 863 825 864 END DO ! ji … … 832 871 !CDIR NODEP 833 872 DO ji = 1, nbrem 834 zji = nind_i(ji)835 zjj = nind_j(ji)836 837 jl1 = zdonor( zji,zjj,jl)873 ii = nind_i(ji) 874 ij = nind_j(ji) 875 876 jl1 = zdonor(ii,ij,jl) 838 877 IF (jl1 .EQ. jl) THEN 839 878 jl2 = jl+1 … … 842 881 ENDIF 843 882 844 zdeice = e_i( zji,zjj,jk,jl1) * zworka(zji,zjj)845 e_i( zji,zjj,jk,jl1) = e_i(zji,zjj,jk,jl1) - zdeice846 e_i( zji,zjj,jk,jl2) = e_i(zji,zjj,jk,jl2) + zdeice883 zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij) 884 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice 885 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 847 886 END DO ! ji 848 887 END DO ! jk … … 860 899 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 861 900 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 862 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl) )) !0 if no ice and 1 if yes901 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 863 902 ELSE 864 903 ht_i(ji,jj,jl) = 0._wp … … 967 1006 zshiftflag = 1 968 1007 zdonor(ji,jj,jl) = jl 969 zdaice(ji,jj,jl) = a_i(ji,jj,jl) 970 zdvice(ji,jj,jl) = v_i(ji,jj,jl) 1008 ! begin TECLIM change 1009 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) 1010 !zdvice(ji,jj,jl) = v_i(ji,jj,jl) 1011 zdaice(ji,jj,jl) = a_i(ji,jj,jl)/2 1012 zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1))/2 1013 ! end TECLIM change 971 1014 ENDIF 972 1015 END DO ! ji -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90
r3625 r4161 26 26 27 27 !!---------------------------------------------------------------------- 28 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)28 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 29 29 !! $Id$ 30 30 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r3791 r4161 41 41 USE agrif_lim2_interp 42 42 #endif 43 #if defined key_bdy 44 USE bdyice_lim 45 #endif 43 46 44 47 IMPLICIT NONE … … 53 56 # include "vectopt_loop_substitute.h90" 54 57 !!---------------------------------------------------------------------- 55 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)58 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 56 59 !! $Id$ 57 60 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 413 416 414 417 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 ) 415 deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 416 !!gm faster to replace the line above with simply: 417 !! deltat(ji,jj) = MAX( delta, creepl ) 418 !!gm end 419 418 ! MV rewriting 419 ! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 420 !!gm faster to replace the line above with simply: 421 !! deltat(ji,jj) = MAX( delta, creepl ) 422 !!gm end 423 deltat(ji,jj) = delta + creepl 424 ! END MV 420 425 !-Calculate stress tensor components zs1 and zs2 421 426 !-at centre of grid cells (see section 3.5 of CICE user's guide). … … 472 477 473 478 CALL lbc_lnk( zs12(:,:), 'F', 1. ) 479 480 !#if defined key_bdy 481 ! ! clem: change zs1, zs2, zs12 at the boundary for each iteration 482 ! CALL bdy_ice_lim_dyn( 2, zs1, zs2, zs12 ) 483 ! CALL lbc_lnk( zs1 (:,:), 'T', 1. ) 484 ! CALL lbc_lnk( zs2 (:,:), 'T', 1. ) 485 ! CALL lbc_lnk( zs12(:,:), 'F', 1. ) 486 !#endif 474 487 475 488 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) … … 520 533 521 534 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 522 #if defined key_agrif 535 #if defined key_agrif && defined key_lim2 523 536 CALL agrif_rhg_lim2( jter, nevp, 'U' ) 524 537 #endif … … 548 561 549 562 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 550 #if defined key_agrif 563 #if defined key_agrif && defined key_lim2 551 564 CALL agrif_rhg_lim2( jter, nevp, 'V' ) 552 565 #endif … … 577 590 578 591 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 579 #if defined key_agrif 592 #if defined key_agrif && defined key_lim2 580 593 CALL agrif_rhg_lim2( jter, nevp , 'V' ) 581 594 #endif … … 608 621 609 622 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 610 #if defined key_agrif 623 #if defined key_agrif && defined key_lim2 611 624 CALL agrif_rhg_lim2( jter, nevp, 'U' ) 612 625 #endif 613 626 614 627 ENDIF 628 629 !#if defined key_bdy 630 ! ! clem: change u_ice and v_ice at the boundary for each iteration 631 ! CALL bdy_ice_lim_dyn( 1 ) 632 !#endif 615 633 616 634 IF(ln_ctl) THEN … … 624 642 ENDIF 625 643 626 ! 644 ! ! ==================== ! 627 645 END DO ! end loop over jter ! 628 646 ! ! ==================== ! 629 630 647 ! 631 648 !------------------------------------------------------------------------------! 632 649 ! 4) Prevent ice velocities when the ice is thin 633 650 !------------------------------------------------------------------------------! 634 ! 635 ! If the ice thickness is below 1cm then ice velocity should equal the 651 !clem : add hminrhg in the namelist 652 ! 653 ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 636 654 ! ocean velocity, 637 655 ! This prevents high velocity when ice is thin … … 641 659 DO ji = fs_2, fs_jpim1 642 660 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 643 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 644 IF ( zdummy .LE. 5.0e-2 ) THEN 661 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 662 zdummy = vt_i(ji,jj) 663 IF ( zdummy .LE. hminrhg ) THEN 645 664 u_ice(ji,jj) = u_oce(ji,jj) 646 665 v_ice(ji,jj) = v_oce(ji,jj) … … 651 670 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 652 671 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 653 #if defined key_agrif 672 #if defined key_agrif && defined key_lim2 654 673 CALL agrif_rhg_lim2( nevp , nevp, 'U' ) 655 674 CALL agrif_rhg_lim2( nevp , nevp, 'V' ) 656 675 #endif 676 #if defined key_bdy 677 ! clem: change u_ice and v_ice at the boundary 678 CALL bdy_ice_lim_dyn( 1 ) 679 #endif 657 680 658 681 DO jj = k_j1+1, k_jpj-1 659 682 DO ji = fs_2, fs_jpim1 660 683 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 661 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 662 IF ( zdummy .LE. 5.0e-2 ) THEN 684 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 685 zdummy = vt_i(ji,jj) 686 IF ( zdummy .LE. hminrhg ) THEN 663 687 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 664 688 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & … … 683 707 !- zds(:,:): shear on northeast corner of grid cells 684 708 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 685 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )686 687 IF ( zdummy .LE. 5.0e-2) THEN709 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 710 zdummy = vt_i(ji,jj) 711 IF ( zdummy .LE. hminrhg ) THEN 688 712 689 713 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & … … 719 743 & - e1v( ji , jj-1 ) * u_ice2(ji ,jj-1) ) / area(ji,jj) 720 744 721 deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) & 722 & + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 723 & ) + creepl 724 745 ! deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) & 746 ! & + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 747 ! & ) + creepl 748 ! MV rewriting 749 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 ) 750 deltat(ji,jj) = delta + creepl 751 ! END MV 752 725 753 ENDIF ! zdummy 726 754 … … 738 766 divu_i (ji,jj) = zdd (ji,jj) 739 767 delta_i(ji,jj) = deltat(ji,jj) 768 ! begin TECLIM change 769 zdst(ji,jj)= ( e2u( ji , jj ) * v_ice1(ji,jj) & 770 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 771 & + e1v( ji , jj ) * u_ice2(ji,jj) & 772 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 740 773 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst(ji,jj) * zdst(ji,jj) ) 774 ! end TECLIM change 741 775 END DO 742 776 END DO 743 CALL lbc_lnk( divu_i (:,:), 'T', 1. ) ! Lateral boundary condition 777 778 ! Lateral boundary condition 779 CALL lbc_lnk( divu_i (:,:), 'T', 1. ) 744 780 CALL lbc_lnk( delta_i(:,:), 'T', 1. ) 781 ! CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 745 782 CALL lbc_lnk( shear_i(:,:), 'T', 1. ) 746 783 -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
r3625 r4161 38 38 39 39 !!---------------------------------------------------------------------- 40 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)40 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 41 41 !! $Id$ 42 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 162 162 CALL iom_rstput( iter, nitrst, numriw, 'v_ice' , v_ice ) 163 163 CALL iom_rstput( iter, nitrst, numriw, 'fsbbq' , fsbbq ) 164 CALL iom_rstput( iter, nitrst, numriw, 'iatte' , iatte ) ! clem modif 165 CALL iom_rstput( iter, nitrst, numriw, 'oatte' , oatte ) ! clem modif 164 166 CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) 165 167 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) … … 340 342 !Control of date 341 343 342 IF( ( nit000 - INT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 ) &344 IF( ( nit000 - NINT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 ) & 343 345 & CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nit000 in ice restart', & 344 346 & ' verify the file or rerun with the value 0 for the', & 345 347 & ' control of time parameter nrstdt' ) 346 IF( INT(zfice) /= nn_fsbc .AND. ABS( nrstdt ) == 1 ) &348 IF( NINT(zfice) /= nn_fsbc .AND. ABS( nrstdt ) == 1 ) & 347 349 & CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nn_fsbc in ice restart', & 348 350 & ' verify the file or rerun with the value 0 for the', & … … 369 371 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 370 372 t_su(:,:,jl) = z2d(:,:) 373 tn_ice (:,:,:) = t_su (:,:,:) 371 374 END DO 372 375 … … 437 440 CALL iom_get( numrir, jpdom_autoglo, 'v_ice' , v_ice ) 438 441 CALL iom_get( numrir, jpdom_autoglo, 'fsbbq' , fsbbq ) 442 CALL iom_get( numrir, jpdom_autoglo, 'iatte' , iatte ) ! clem modif 443 CALL iom_get( numrir, jpdom_autoglo, 'oatte' , oatte ) ! clem modif 439 444 CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i ) 440 445 CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i ) … … 563 568 END DO 564 569 ! 565 CALL iom_close( numrir )570 !clem CALL iom_close( numrir ) 566 571 ! 567 572 CALL wrk_dealloc( nlay_i, zs_zero ) -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r4148 r4161 10 10 !! ! + simplification of the ice-ocean stress calculation 11 11 !! 3.4 ! 2011-02 (G. Madec) dynamical allocation 12 !! - ! 2012 (D. Iovino) salt flux change 13 !! - ! 2012-05 (C. Rousset) add penetration solar flux 12 14 !! 3.5 ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass 13 15 !!---------------------------------------------------------------------- … … 35 37 USE prtctl ! Print control 36 38 USE cpl_oasis3, ONLY : lk_cpl 39 USE traqsr ! clem: add penetration of solar flux into the calculation of heat budget 37 40 USE oce, ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n 38 41 USE dom_ice, ONLY : tms … … 57 60 # include "vectopt_loop_substitute.h90" 58 61 !!---------------------------------------------------------------------- 59 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)62 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 60 63 !! $Id$ 61 64 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 99 102 INTEGER, INTENT(in) :: kt ! number of iteration 100 103 ! 101 INTEGER :: ji, jj ! dummy loop indices104 INTEGER :: ji, jj, jl ! dummy loop indices 102 105 INTEGER :: ierr, ifvt, i1mfr, idfr ! local integer 103 106 INTEGER :: iflt, ial , iadv , ifral, ifrdv ! - - … … 106 109 REAL(wp) :: zfcm1 , zfcm2 ! - - 107 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb, zalbp ! 2D/3D workspace 111 REAL(wp) :: zzfcm1, zfscmbq ! clem: for light penetration 108 112 !!--------------------------------------------------------------------- 109 113 … … 119 123 DO ji = 1, jpi 120 124 zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 121 ifvt = zinda * MAX( rzero , SIGN( rone, - phicif(ji,jj) ) ) !subscripts are bad here122 i1mfr = 1.0 - MAX( rzero , SIGN( rone , - ( at_i(ji,jj)) ) )125 ifvt = zinda * MAX( rzero , SIGN( rone, - phicif(ji,jj) ) ) !subscripts are bad here 126 i1mfr = 1.0 - MAX( rzero , SIGN( rone , - at_i(ji,jj) ) ) 123 127 idfr = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) ) 124 128 iflt = zinda * (1 - i1mfr) * (1 - ifvt ) … … 141 145 142 146 ! computation the solar flux at ocean surface 143 zfcm1 = pfrld(ji,jj) * qsr(ji,jj) + ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) 147 IF (lk_cpl) THEN ! be carfeful: not being tested yet 148 ! original line 149 !zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) 150 ! new line to include solar penetration (not tested) 151 zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 152 DO jl = 1, jpl 153 zfcm1 = zfcm1 - qsr_ice(ji,jj,jl) * a_i(ji,jj,jl) 154 END DO 155 ELSE 156 zfcm1 = pfrld(ji,jj) * qsr(ji,jj) + & 157 & ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 158 ENDIF 144 159 ! fstric Solar flux transmitted trough the ice 145 160 ! qsr Net short wave heat flux on free ocean 146 161 ! new line 147 fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) 162 fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 163 164 ! solar flux and fscmbq with light penetration (clem) 165 zzfcm1 = pfrld(ji,jj) * qsr(ji,jj) * oatte(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) 166 zfscmbq = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) 148 167 149 168 ! computation the non solar heat flux at ocean surface 150 zfcm2 = - z fcm1 & ! ???151 & + iflt * fscmbq(ji,jj)& ! total ablation: heat given to the ocean169 zfcm2 = - zzfcm1 & ! 170 & + iflt * zfscmbq & ! total ablation: heat given to the ocean 152 171 & + ifral * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice & 153 172 & + ifrdv * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * r1_rdtice & … … 170 189 ! ! fdtcn : turbulent oceanic heat flux 171 190 172 !!gm this IF prevents the vertorisation of the whole loop173 IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN174 WRITE(numout,*) ' lim_sbc : heat fluxes '175 WRITE(numout,*) ' qsr : ', qsr(jiindx,jjindx)176 WRITE(numout,*) ' pfrld : ', pfrld(jiindx,jjindx)177 WRITE(numout,*) ' fstric : ', fstric (jiindx,jjindx)178 WRITE(numout,*)179 WRITE(numout,*) ' qns : ', qns(jiindx,jjindx)180 WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)181 WRITE(numout,*) ' ifral : ', ifral182 WRITE(numout,*) ' ial : ', ial183 WRITE(numout,*) ' qcmif : ', qcmif(jiindx,jjindx)184 WRITE(numout,*) ' qldif : ', qldif(jiindx,jjindx)185 186 187 WRITE(numout,*) ' ifrdv : ', ifrdv188 WRITE(numout,*) ' qfvbq : ', qfvbq(jiindx,jjindx)189 WRITE(numout,*) ' qdtcn : ', qdtcn(jiindx,jjindx)190 191 192 WRITE(numout,*) ' '193 WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)194 WRITE(numout,*) ' fhmec : ', fhmec(jiindx,jjindx)195 WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx)196 WRITE(numout,*) ' fhbri : ', fhbri(jiindx,jjindx)197 WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)198 ENDIF199 !!gm end191 !!gm this IF prevents the vertorisation of the whole loop 192 ! IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 193 ! WRITE(numout,*) ' lim_sbc : heat fluxes ' 194 ! WRITE(numout,*) ' qsr : ', qsr(jiindx,jjindx) 195 ! WRITE(numout,*) ' pfrld : ', pfrld(jiindx,jjindx) 196 ! WRITE(numout,*) ' fstric : ', fstric (jiindx,jjindx) 197 ! WRITE(numout,*) 198 ! WRITE(numout,*) ' qns : ', qns(jiindx,jjindx) 199 ! WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx) 200 ! WRITE(numout,*) ' ifral : ', ifral 201 ! WRITE(numout,*) ' ial : ', ial 202 ! WRITE(numout,*) ' qcmif : ', qcmif(jiindx,jjindx) 203 ! WRITE(numout,*) ' qldif : ', qldif(jiindx,jjindx) 204 ! !WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice 205 ! !WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice 206 ! WRITE(numout,*) ' ifrdv : ', ifrdv 207 ! WRITE(numout,*) ' qfvbq : ', qfvbq(jiindx,jjindx) 208 ! WRITE(numout,*) ' qdtcn : ', qdtcn(jiindx,jjindx) 209 ! !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice 210 ! !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice 211 ! WRITE(numout,*) ' ' 212 ! WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx) 213 ! WRITE(numout,*) ' fhmec : ', fhmec(jiindx,jjindx) 214 ! WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx) 215 ! WRITE(numout,*) ' fhbri : ', fhbri(jiindx,jjindx) 216 ! WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 217 ! ENDIF 218 !!gm end 200 219 END DO 201 220 END DO … … 218 237 219 238 ! computing freshwater exchanges at the ice/ocean interface 220 zemp = emp(ji,jj) * ( 1.0 - at_i(ji,jj) ) & ! evaporation over oceanic fraction 221 & - tprecip(ji,jj) * at_i(ji,jj) & ! all precipitation reach the ocean 222 & + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) ) & ! except solid precip intercepted by sea-ice 223 & - fmmec(ji,jj) ! snow falling when ridging 239 IF (lk_cpl) THEN 240 zemp = - emp_tot(ji,jj) + emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) ) & ! 241 & - rdm_snw(ji,jj) / rdt_ice 242 ELSE 243 zemp = emp(ji,jj) * ( 1.0 - at_i(ji,jj) ) & ! evaporation over oceanic fraction 244 & - tprecip(ji,jj) * at_i(ji,jj) & ! all precipitation reach the ocean 245 & + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) ) & ! except solid precip intercepted by sea-ice 246 & - fmmec(ji,jj) ! snow falling when ridging 247 ENDIF 224 248 225 249 ! mass flux at the ocean/ice interface (sea ice fraction) … … 370 394 !! ** input : Namelist namicedia 371 395 !!------------------------------------------------------------------- 396 REAL(wp) :: zsum, zarea 372 397 ! 373 398 INTEGER :: ji, jj ! dummy loop indices … … 390 415 END WHERE 391 416 ENDIF 417 ! clem modif 418 iatte(:,:) = 1._wp 419 oatte(:,:) = 1._wp 420 ! 392 421 ! ! embedded sea ice 393 422 IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass … … 435 464 ENDIF 436 465 ! 466 !!? IF( .NOT. ln_rstart ) THEN ! delete the initial ssh below sea-ice area 467 !!? ! 468 !!? zarea = glob_sum( e1e2t(:,:) ) ! interior global domain surface 469 !!? zsum = glob_sum( e1e2t(:,:) * ( snwice_mass(:,:) ) ) / zarea * r1_rau0 470 !!? sshn(:,:) = sshn(:,:) - zsum 471 !!? sshb(:,:) = sshb(:,:) - zsum 472 !!? ENDIF 473 ! 474 437 475 END SUBROUTINE lim_sbc_init 438 476 -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90
r3625 r4161 20 20 21 21 !!---------------------------------------------------------------------- 22 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2010)22 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 23 23 !! $Id$ 24 24 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4147 r4161 11 11 !! 3.3 ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas) 12 12 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 13 !! - ! 2012-05 (C. Rousset) add penetration solar flux 13 14 !!---------------------------------------------------------------------- 14 15 #if defined key_lim3 … … 40 41 USE prtctl ! Print control 41 42 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 43 USE timing ! Timing 42 44 43 45 IMPLICIT NONE … … 92 94 REAL(wp) :: zfntlat, zpareff, zareamin, zcoef ! - - 93 95 REAL(wp), POINTER, DIMENSION(:,:) :: zqlbsbq ! link with lead energy budget qldif 96 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 97 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 94 98 !!------------------------------------------------------------------- 99 IF( nn_timing == 1 ) CALL timing_start('limthd') 95 100 96 101 CALL wrk_alloc( jpi, jpj, zqlbsbq ) 97 102 103 ! ------------------------------- 104 !- check conservation (C Rousset) 105 IF (ln_limdiahsb) THEN 106 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 107 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 108 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 109 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 110 ENDIF 111 !- check conservation (C Rousset) 112 ! ------------------------------- 113 98 114 !------------------------------------------------------------------------------! 99 115 ! 1) Initialization of diagnostic variables ! … … 109 125 DO ji = 1, jpi 110 126 !Energy of melting q(S,T) [J.m-3] 111 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * nlay_i127 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * REAL( nlay_i ) 112 128 !0 if no ice and 1 if yes 113 129 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) ) ) … … 121 137 DO ji = 1, jpi 122 138 !Energy of melting q(S,T) [J.m-3] 123 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * nlay_s139 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL( nlay_s ) 124 140 !0 if no ice and 1 if yes 125 141 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) ) ) … … 134 150 ! 1.3) Set some dummies to 0 135 151 !----------------------------- 136 rdvosif(:,:) = 0.e0 ! variation of ice volume at surface137 rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom138 fdvolif(:,:) = 0.e0 ! total variation of ice volume139 rdvonif(:,:) = 0.e0 ! lateral variation of ice volume140 fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice141 ffltbif(:,:) = 0.e0 ! linked with fstric142 qfvbq (:,:) = 0.e0 ! linked with fstric143 rdm_snw(:,:) = 0.e0 ! variation of snow mass per unit area144 rdm_ice(:,:) = 0.e0 ! variation of ice mass per unit area145 hicifp (:,:) = 0.e0 ! daily thermodynamic ice production.146 sfx_bri(:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean147 fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean148 sfx_thd(:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay152 !clem rdvosif(:,:) = 0.e0 ! variation of ice volume at surface 153 !clem rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom 154 !clem fdvolif(:,:) = 0.e0 ! total variation of ice volume 155 !clem rdvonif(:,:) = 0.e0 ! lateral variation of ice volume 156 !clem fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice 157 !clem ffltbif(:,:) = 0.e0 ! linked with fstric 158 !clem qfvbq (:,:) = 0.e0 ! linked with fstric 159 !clem rdm_snw(:,:) = 0.e0 ! variation of snow mass per unit area 160 !clem rdm_ice(:,:) = 0.e0 ! variation of ice mass per unit area 161 !clem hicifp (:,:) = 0.e0 ! daily thermodynamic ice production. 162 !clem sfx_bri(:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean 163 !clem fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean 164 !clem sfx_thd(:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay 149 165 150 166 !----------------------------------- … … 165 181 !CDIR NOVERRCHK 166 182 DO ji = 1, jpi 167 zthsnice = SUM( ht_s(ji,jj,1:jpl) ) + SUM( ht_i(ji,jj,1:jpl) )168 zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) )169 183 phicif(ji,jj) = vt_i(ji,jj) 170 184 pfrld(ji,jj) = 1.0 - at_i(ji,jj) 171 zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )185 zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) ) 172 186 ! 173 187 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 180 194 181 195 ! here the drag will depend on ice thickness and type (0.006) 182 fdtcn(ji,jj) = zind b* rau0 * rcp * 0.006 * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) )196 fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) ) 183 197 ! also category dependent 184 198 ! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead 185 qdtcn(ji,jj) = zind b* fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice199 qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 186 200 ! 187 201 ! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) 188 202 ! ! caution: exponent betas used as more snow can fallinto leads 189 203 qldif(ji,jj) = tms(ji,jj) * rdt_ice * ( & 190 & pfrld(ji,jj) * ( qsr(ji,jj) & ! solar heat204 & pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif 191 205 & + qns(ji,jj) & ! non solar heat 192 206 & + fdtcn(ji,jj) & ! turbulent ice-ocean heat 193 & + fsbbq(ji,jj) * ( 1.0 - zind b) ) & ! residual heat from previous step207 & + fsbbq(ji,jj) * ( 1.0 - zinda ) ) & ! residual heat from previous step 194 208 & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus ) ! latent heat of sprecip melting 195 209 ! … … 206 220 ! 207 221 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 208 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda )222 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) 209 223 ! 210 224 ! oceanic heat flux (limthd_dh) 211 fbif (ji,jj) = zind b* ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) )225 fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 212 226 ! 213 227 END DO … … 294 308 CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb), qfvbq , jpi, jpj, npb(1:nbpb) ) 295 309 310 CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) ) ! clem modif 311 CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) ) ! clem modif 296 312 !-------------------------------- 297 313 ! 4.3) Thermodynamic processes … … 411 427 ! 5.4) Diagnostic thermodynamic growth rates 412 428 !-------------------------------------------- 413 d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes414 429 !clem@useless d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes 430 !clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 415 431 416 432 IF( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) … … 448 464 ENDIF 449 465 ! 466 ! ------------------------------- 467 !- check conservation (C Rousset) 468 IF (ln_limdiahsb) THEN 469 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 470 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 471 472 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 473 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 474 475 zchk_vmin = glob_min(v_i) 476 zchk_amax = glob_max(SUM(a_i,dim=3)) 477 zchk_amin = glob_min(a_i) 478 479 IF(lwp) THEN 480 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limthd) = ',(zchk_v_i * rday) 481 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * rday) 482 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limthd) = ',(zchk_vmin * 1.e-3) 483 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limthd) = ',zchk_amax 484 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limthd) = ',zchk_amin 485 ENDIF 486 ENDIF 487 !- check conservation (C Rousset) 488 ! ------------------------------- 489 ! 450 490 CALL wrk_dealloc( jpi, jpj, zqlbsbq ) 451 491 ! 492 IF( nn_timing == 1 ) CALL timing_stop('limthd') 452 493 END SUBROUTINE lim_thd 453 494 … … 472 513 DO jk = 1, nlay_i ! total q over all layers, ice [J.m-2] 473 514 DO ji = kideb, kiut 474 etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i515 etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 475 516 eti (ji,jl) = eti(ji,jl) + etilayer(ji,jk) 476 517 END DO 477 518 END DO 478 519 DO ji = kideb, kiut ! total q over all layers, snow [J.m-2] 479 ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / nlay_s520 ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / REAL( nlay_s ) 480 521 END DO 481 522 ! … … 498 539 499 540 INTEGER :: ji, jk ! loop indices 500 INTEGER :: zji, zjj541 INTEGER :: ii, ij 501 542 INTEGER :: numce ! number of points for which conservation is violated 502 543 REAL(wp) :: meance ! mean conservation error … … 521 562 !---------------------------------------- 522 563 DO ji = kideb, kiut 523 zji = MOD( npb(ji) - 1 , jpi ) + 1524 zjj = ( npb(ji) - 1 ) / jpi + 1564 ii = MOD( npb(ji) - 1 , jpi ) + 1 565 ij = ( npb(ji) - 1 ) / jpi + 1 525 566 fatm (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) 526 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc( zji,zjj,jl)567 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(ii,ij,jl) 527 568 END DO 528 569 … … 579 620 IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. & 580 621 ( cons_error(ji,jl) .GT. max_cons_err ) ) THEN 581 zji = MOD( npb(ji) - 1, jpi ) + 1582 zjj = ( npb(ji) - 1 ) / jpi + 1622 ii = MOD( npb(ji) - 1, jpi ) + 1 623 ij = ( npb(ji) - 1 ) / jpi + 1 583 624 ! 584 625 WRITE(numout,*) ' alerte 1 ' … … 586 627 WRITE(numout,*) ' heat diffusion in the ice ' 587 628 WRITE(numout,*) ' Category : ', jl 588 WRITE(numout,*) ' zji , zjj : ', zji, zjj589 WRITE(numout,*) ' lat, lon : ', gphit( zji,zjj), glamt(zji,zjj)629 WRITE(numout,*) ' ii , ij : ', ii, ij 630 WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij) 590 631 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 591 632 WRITE(numout,*) ' surf_error : ', surf_error(ji,jl) … … 615 656 WRITE(numout,*) ' fc_bo : ', - fc_bo_i (ji) 616 657 WRITE(numout,*) ' foc : ', fbif_1d(ji) 617 WRITE(numout,*) ' fstroc : ', fstroc ( zji,zjj,jl)658 WRITE(numout,*) ' fstroc : ', fstroc (ii,ij,jl) 618 659 WRITE(numout,*) ' i0 : ', i0(ji) 619 660 WRITE(numout,*) ' qsr_ice : ', (1.0-i0(ji))*qsr_ice_1d(ji) … … 651 692 ! 652 693 INTEGER :: ji ! loop indices 653 INTEGER :: zji, zjj, numce ! local integers694 INTEGER :: ii, ij, numce ! local integers 654 695 REAL(wp) :: meance, max_cons_err !local scalar 655 696 !!--------------------------------------------------------------------- … … 669 710 !---------------------------------------- 670 711 DO ji = kideb, kiut 671 zji = MOD( npb(ji) - 1 , jpi ) + 1672 zjj = ( npb(ji) - 1 ) / jpi + 1712 ii = MOD( npb(ji) - 1 , jpi ) + 1 713 ij = ( npb(ji) - 1 ) / jpi + 1 673 714 674 715 fatm (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji) ! total heat flux 675 sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc( zji,zjj,jl)716 sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(ii,ij,jl) 676 717 cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 677 718 END DO … … 706 747 DO ji = kideb, kiut 707 748 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 708 zji = MOD( npb(ji) - 1, jpi ) + 1709 zjj = ( npb(ji) - 1 ) / jpi + 1749 ii = MOD( npb(ji) - 1, jpi ) + 1 750 ij = ( npb(ji) - 1 ) / jpi + 1 710 751 ! 711 752 WRITE(numout,*) ' alerte 1 - category : ', jl 712 753 WRITE(numout,*) ' Untolerated conservation error after limthd_ent ' 713 WRITE(numout,*) ' zji , zjj : ', zji, zjj714 WRITE(numout,*) ' lat, lon : ', gphit( zji,zjj), glamt(zji,zjj)754 WRITE(numout,*) ' ii , ij : ', ii, ij 755 WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij) 715 756 WRITE(numout,*) ' * ' 716 757 WRITE(numout,*) ' Ftotal : ', sum_fluxq(ji,jl) … … 724 765 WRITE(numout,*) ' foce : ', fbif_1d(ji) 725 766 WRITE(numout,*) ' fres : ', ftotal_fin(ji) 726 WRITE(numout,*) ' fhbri : ', fhbricat( zji,zjj,jl)767 WRITE(numout,*) ' fhbri : ', fhbricat(ii,ij,jl) 727 768 WRITE(numout,*) ' * ' 728 769 WRITE(numout,*) ' Heat contents --- : ' … … 793 834 INTEGER :: ios ! Local integer output status for namelist read 794 835 NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, & 795 & hicmin, hiclim, amax ,&836 & hicmin, hiclim, & 796 837 & sbeta , parlat, hakspl, hibspl, exld, & 797 838 & hakdif, hnzst , thth , parsub, alphs, betas, & … … 825 866 WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin 826 867 WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim 827 WRITE(numout,*)' maximum lead fraction amax = ', amax828 868 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 829 869 WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r3808 r4161 6 6 !! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D 7 7 !! ! 2005-06 (M. Vancoppenolle) 3D version 8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm snif & rdmicif8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw & rdm_ice 9 9 !! 3.4 ! 2011-02 (G. Madec) dynamical allocation 10 10 !! 3.5 ! 2012-10 (G. Madec & co) salt flux + bug fixes … … 39 39 40 40 !!---------------------------------------------------------------------- 41 !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)41 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 42 42 !! $Id$ 43 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 73 73 !! 74 74 INTEGER :: ji , jk ! dummy loop indices 75 INTEGER :: ii, ij ! 2D corresponding indices to ji75 INTEGER :: ii, ij ! 2D corresponding indices to ji 76 76 INTEGER :: isnow ! switch for presence (1) or absence (0) of snow 77 77 INTEGER :: isnowic ! snow ice formation not … … 84 84 REAL(wp) :: zdhnm, zhnnew, zhisn, zihic, zzc ! 85 85 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 86 REAL(wp) :: zds ! increment of bottom ice salinity87 86 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads 88 87 REAL(wp) :: zsm_snowice ! snow-ice salinity … … 107 106 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_pre ! snow precipitation 108 107 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_sub ! snow sublimation 109 REAL(wp), POINTER, DIMENSION(:) :: zsfx_melt ! salt flux due to ice melt110 108 111 109 REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah … … 120 118 REAL(wp), POINTER, DIMENSION(:,:) :: zqt_i_lay ! total ice heat content 121 119 120 ! mass and salt flux (clem) 121 REAL(wp) :: zdvres, zdvsur, zdvbot 122 REAL(wp), POINTER, DIMENSION(:) :: zviold, zvsold ! old ice volume... 123 122 124 ! Heat conservation 123 125 INTEGER :: num_iter_max, numce_dh … … 128 130 129 131 CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 130 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, z sfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )132 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 131 133 CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i ) 132 134 CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) 133 135 134 zsfx_melt (:) = 0._wp 136 CALL wrk_alloc( jpij, zviold, zvsold ) ! clem 137 135 138 ftotal_fin(:) = 0._wp 136 139 zfdt_init (:) = 0._wp 137 140 zfdt_final(:) = 0._wp 138 141 142 dh_i_surf (:) = 0._wp 143 dh_i_bott (:) = 0._wp 144 dh_snowice(:) = 0._wp 145 139 146 DO ji = kideb, kiut 140 147 old_ht_i_b(ji) = ht_i_b(ji) 141 148 old_ht_s_b(ji) = ht_s_b(ji) 149 zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem 150 zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem 142 151 END DO 143 152 ! … … 164 173 ! 165 174 DO ji = kideb, kiut ! Layer thickness 166 zh_i(ji) = ht_i_b(ji) / nlay_i167 zh_s(ji) = ht_s_b(ji) / nlay_s175 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 176 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 168 177 END DO 169 178 ! … … 171 180 DO jk = 1, nlay_s 172 181 DO ji = kideb, kiut 173 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s182 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) 174 183 END DO 175 184 END DO … … 178 187 DO jk = 1, nlay_i 179 188 DO ji = kideb, kiut 180 zzc = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i189 zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 181 190 zqt_i(ji) = zqt_i(ji) + zzc 182 191 zqt_i_lay(ji,jk) = zzc … … 244 253 zhn = 1.0 - MAX( zzero , SIGN( zone , - zhsnew ) ) 245 254 ht_s_b(ji) = MAX( zzero , zhsnew ) 255 ! we recompute dh_s_tot (clem) 256 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 246 257 ! Volume and mass variations of snow 247 258 dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) 248 259 dvsbq_1d (ji) = MIN( zzero, dvsbq_1d(ji) ) 249 rdm_snw_1d(ji) = rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji)260 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 250 261 END DO ! ji 251 262 … … 254 265 !-------------------------- 255 266 DO ji = kideb, kiut 256 dh_i_surf(ji) = 0._wp257 267 z_f_surf (ji) = zqfont_su(ji) * r1_rdtice ! heat conservation test 258 268 zdq_i (ji) = 0._wp … … 272 282 zdq_i (ji ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 273 283 ! 274 ! ! contribution to ice-ocean salt flux 275 zsfx_melt(ji) = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 284 ! clem 285 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 286 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 276 287 END DO 277 288 END DO … … 334 345 DO ji = kideb,kiut 335 346 q_s_b (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 336 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s! heat conservation347 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) ! heat conservation 337 348 END DO 338 349 END DO … … 375 386 ! Basal growth rate = - F*dt / q 376 387 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 388 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 377 389 ENDIF 378 390 END DO … … 416 428 zfracs = zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & 417 429 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 418 z ds = zfracs * sss_m(ii,ij) - s_i_new(ji)430 zfracs = MIN( 0.5 , zfracs ) 419 431 s_i_new(ji) = zfracs * sss_m(ii,ij) 420 432 ENDIF ! fc_bo_i … … 425 437 DO ji = kideb, kiut 426 438 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN 427 ! New ice salinity must not exceed 15psu439 ! New ice salinity must not exceed 20 psu 428 440 s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 429 441 ! Metling point in K … … 437 449 ! Salinity update 438 450 ! entrapment during bottom growth 439 dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) ) & 440 & / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji) 451 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 441 452 ENDIF ! heat budget 442 453 END DO … … 476 487 zdq_i (ji ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 477 488 ENDIF 478 ! contribution to salt flux 479 zsfx_melt(ji) = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 489 ! clem: contribution to salt flux 490 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 491 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 480 492 ENDIF 481 493 END DO ! ji … … 528 540 ELSE ; zdhbf = dh_i_bott(ji) 529 541 ENDIF 542 zdvres = zdhbf - dh_i_bott(ji) 543 dh_i_bott(ji) = zdhbf 544 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 530 545 ! ! excessive energy is sent to lateral ablation 531 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) & 532 & * ( zdhbf - dh_i_bott(ji) ) * r1_rdtice 533 dh_i_bott(ji) = zdhbf 534 ! !since ice volume is only used for outputs, we keep it global for all categories 535 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 536 ! !new ice thickness 537 zhgnew (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 538 ! ! diagnostic ( bottom ice growth ) 539 ii = MOD( npb(ji) - 1, jpi ) + 1 540 ij = ( npb(ji) - 1 ) / jpi + 1 541 diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 542 diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 543 diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 546 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres * r1_rdtice 544 547 END DO 545 548 … … 552 555 ! Adapt the remaining energy if too much ice melts 553 556 !-------------------------------------------------- 554 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) ! =1 if ice 555 ! 0 if no more ice 556 zhgnew (ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 557 ! remaining heat 557 zdvres = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) ) 558 zdvsur = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first 559 zdvbot = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom 560 dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem 561 dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem 562 563 ! new ice thickness (clem) 564 zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 565 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 566 zhgnew(ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 567 568 ! !since ice volume is only used for outputs, we keep it global for all categories 569 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 570 571 ! remaining heat 558 572 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) 559 573 … … 569 583 ht_s_b(ji) = MAX( zzero , zhnfi ) 570 584 zqt_s(ji) = zqt_s(ji) * ht_s_b(ji) 585 ! we recompute dh_s_tot (clem) 586 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 571 587 572 588 ! Mass variations of ice and snow … … 579 595 ! 580 596 ! ! mass variation cumulated over category 581 rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s ! snow582 rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i ! ice597 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s ! snow 598 !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i ! ice 583 599 584 600 ! Remaining heat to the ocean … … 586 602 focea(ji) = - zfdt_final(ji) * r1_rdtice ! focea is in W.m-2 * dt 587 603 604 ! residual salt flux (clem) 605 !-------------------------- 606 ! surface 607 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvsur * rhoic * r1_rdtice 608 ! bottom 609 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting 610 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 611 ELSE ! growth 612 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 613 ENDIF 614 ! 615 ! diagnostic ( bottom ice growth ) 616 ii = MOD( npb(ji) - 1, jpi ) + 1 617 ij = ( npb(ji) - 1 ) / jpi + 1 618 diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 619 diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 620 diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 588 621 END DO 589 622 … … 591 624 592 625 !--------------------------- 593 ! Salt flux andheat fluxes626 ! heat fluxes 594 627 !--------------------------- 595 628 DO ji = kideb, kiut 596 629 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) ! =1 if ice 597 !598 ! Salt flux599 sfx_thd_1d(ji) = sfx_thd_1d(ji) + zihgnew * zsfx_melt(ji) &600 & - (1.0 - zihgnew) * zfmass_i (ji) * sm_i_b(ji) * r1_rdtice601 630 ! 602 631 ! Heat flux … … 646 675 dmgwi_1d (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 647 676 648 ! All snow is thrown in the ocean, and seawater is taken to replace the volume 649 rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic * ( 1. - rhosn / rhoic ) 650 rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 677 !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 678 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 651 679 652 680 ! Equivalent salt flux (1) Snow-ice formation component … … 658 686 ELSE ; zsm_snowice = sm_i_b(ji) 659 687 ENDIF 660 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice661 !662 688 ! entrapment during snow ice formation 663 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 664 isnowic = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch 665 IF( num_sal == 2 ) & 666 dsm_i_si_1d(ji) = ( zsm_snowice * dh_snowice(ji) & 667 & + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 ) & 668 & - sm_i_b(ji) ) * isnowic 689 ! clem: new salinity difference stored (to be used in limthd_ent.F90) 690 IF ( num_sal == 2 ) THEN 691 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) ) 692 ! salinity dif due to snow-ice formation 693 dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 694 ! salinity dif due to bottom growth 695 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0._wp ) THEN 696 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 697 ENDIF 698 ENDIF 669 699 670 700 ! Actualize new snow and ice thickness. … … 680 710 diag_sni_gr(ii,ij) = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 681 711 ! 712 ! salt flux 713 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 714 !-------------------------------- 715 ! Update mass fluxes (clem) 716 !-------------------------------- 717 rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic 718 rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn 719 682 720 END DO !ji 683 721 ! 684 722 CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 685 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, z sfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )723 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 686 724 CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) 687 725 CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) 726 ! 727 CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem 688 728 ! 689 729 END SUBROUTINE lim_thd_dh -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r3808 r4161 10 10 !! ! 04-2007 (M. Vancoppenolle) Energy conservation 11 11 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 12 !! - ! 2012-05 (C. Rousset) add penetration solar flux 12 13 !!---------------------------------------------------------------------- 13 14 #if defined key_lim3 … … 34 35 35 36 !!---------------------------------------------------------------------- 36 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)37 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 37 38 !! $Id$ 38 39 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 156 157 DO ji = kideb , kiut 157 158 ! is there snow or not 158 isnow(ji)= INT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) )) )159 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) ) 159 160 ! surface temperature of fusion 160 161 !!gm ??? ztfs(ji) = rtt !!!???? 161 ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt162 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 162 163 ! layer thickness 163 zh_i(ji) = ht_i_b(ji) / nlay_i164 zh_s(ji) = ht_s_b(ji) / nlay_s164 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 165 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 165 166 END DO 166 167 … … 174 175 DO layer = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 175 176 DO ji = kideb , kiut 176 z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s177 z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s ) 177 178 END DO 178 179 END DO … … 180 181 DO layer = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 181 182 DO ji = kideb , kiut 182 z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i183 z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i ) 183 184 END DO 184 185 END DO … … 201 202 DO ji = kideb , kiut 202 203 ! switches 203 isnow(ji) = INT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) )) )204 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ) ) 204 205 ! hs > 0, isnow = 1 205 206 zhsu (ji) = hnzst ! threshold for the computation of i0 206 207 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) ) 207 208 208 i0(ji) = ( 1._wp- isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )209 i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 209 210 !fr1_i0_1d = i0 for a thin ice surface 210 211 !fr1_i0_2d = i0 for a thick ice surface … … 243 244 244 245 DO ji = kideb, kiut ! ice initialization 245 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp- isnow(ji) )246 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1 - isnow(ji) ) 246 247 END DO 247 248 … … 256 257 257 258 DO ji = kideb, kiut ! Radiation transmitted below the ice 258 fstbif_1d(ji) = fstbif_1d(ji) + zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji)259 fstbif_1d(ji) = fstbif_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif 259 260 END DO 260 261 … … 264 265 ii = MOD( npb(ji) - 1 , jpi ) + 1 265 266 ij = ( npb(ji) - 1 ) / jpi + 1 266 fstroc(ii,ij,jl) = zradtr_i(ji,nlay_i)267 fstroc(ii,ij,jl) = iatte_1d(ji) * zradtr_i(ji,nlay_i) ! clem modif 267 268 END DO 268 269 ! +++++ … … 376 377 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & 377 378 (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 378 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)* isnow(ji) &379 + zkappa_i(ji,0)* (1.0-isnow(ji))379 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 380 + zkappa_i(ji,0)*REAL( 1 - isnow(ji) ) 380 381 END DO 381 382 ! … … 658 659 t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 659 660 * t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) & 660 * MAX(0.0,SIGN(1.0,ht_s_b(ji) -zeps))661 * MAX(0.0,SIGN(1.0,ht_s_b(ji))) 661 662 662 663 ! surface temperature 663 isnow(ji) = INT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) ) ) )664 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) ) ) ) 664 665 ztsuoldit(ji) = t_su_b(ji) 665 IF( t_su_b(ji) < ztfs(ji) ) 666 t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( isnow(ji)*t_s_b(ji,1) &667 & + (1.0-isnow(ji))*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))666 IF( t_su_b(ji) < ztfs(ji) ) & 667 t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) & 668 & + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 668 669 END DO 669 670 ! … … 721 722 #endif 722 723 ! ! surface ice conduction flux 723 isnow(ji) = INT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) )724 fc_su(ji) = - isnow(ji)* zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) &725 & - ( 1._wp- isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji))724 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) ) 725 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) & 726 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji)) 726 727 ! ! bottom ice conduction flux 727 728 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) … … 734 735 DO ji = kideb, kiut 735 736 ! Upper snow value 736 fc_s(ji,0) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) )737 fc_s(ji,0) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) ) 737 738 ! Bott. snow value 738 fc_s(ji,1) = - isnow(ji)* zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) )739 fc_s(ji,1) = - REAL( isnow(ji) ) * zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) ) 739 740 END DO 740 741 DO ji = kideb, kiut ! Upper ice layer 741 fc_i(ji,0) = - isnow(ji) * & ! interface flux if there is snow742 fc_i(ji,0) = - REAL( isnow(ji) ) * & ! interface flux if there is snow 742 743 ( zkappa_i(ji,0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) & 743 - ( 1.0- isnow(ji) ) * ( zkappa_i(ji,0) * &744 - REAL( 1 - isnow(ji) ) * ( zkappa_i(ji,0) * & 744 745 zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not 745 746 END DO -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r3625 r4161 44 44 45 45 !!---------------------------------------------------------------------- 46 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)46 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 47 47 !! $Id$ 48 48 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 75 75 76 76 INTEGER :: ji,jk ! dummy loop indices 77 INTEGER :: zji, zjj , & ! dummy indices77 INTEGER :: ii, ij , & ! dummy indices 78 78 ntop0 , & ! old layer top index 79 79 nbot1 , & ! new layer bottom index … … 145 145 146 146 DO ji = kideb, kiut 147 zh_i(ji) = old_ht_i_b(ji) / nlay_i148 zh_s(ji) = old_ht_s_b(ji) / nlay_s147 zh_i(ji) = old_ht_i_b(ji) / REAL( nlay_i ) 148 zh_s(ji) = old_ht_s_b(ji) / REAL( nlay_s ) 149 149 END DO 150 150 … … 166 166 DO jk = 1, nlays0 167 167 DO ji = kideb, kiut 168 snind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))) &169 + snind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))))168 snind(ji) = jk * NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)))) & 169 + snind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji))))) 170 170 zdeltah(ji)= zdeltah(ji) + zh_s(ji) 171 171 END DO ! ji … … 175 175 ! 0 if not 176 176 DO ji = kideb, kiut 177 snswi(ji) = MAX(0, INT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji)))))177 snswi(ji) = MAX(0,NINT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji))))) 178 178 END DO ! ji 179 179 … … 190 190 DO jk = 1, nlayi0 191 191 DO ji = kideb, kiut 192 icsuind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))) &193 + icsuind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))))192 icsuind(ji) = jk * NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)))) & 193 + icsuind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji))))) 194 194 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 195 195 END DO ! ji … … 200 200 ! 0 if not 201 201 DO ji = kideb, kiut 202 icsuswi(ji) = MAX(0, INT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) )202 icsuswi(ji) = MAX(0,NINT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) ) 203 203 ENDDO 204 204 … … 216 216 DO jk = nlayi0, 1, -1 217 217 DO ji = kideb, kiut 218 icboind(ji) = (nlayi0+1-jk) * INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))) &219 & + icboind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))))218 icboind(ji) = (nlayi0+1-jk) * NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)))) & 219 & + icboind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji))))) 220 220 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 221 221 END DO … … 232 232 ! 0 if ablation is on the way 233 233 DO ji = kideb, kiut 234 icboswi(ji) = MAX(0, INT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji)))))234 icboswi(ji) = MAX(0,NINT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji))))) 235 235 END DO 236 236 … … 248 248 DO ji = kideb, kiut 249 249 snicind(ji) = (nlays0+1-jk) & 250 * INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))) + snicind(ji) &251 * (1 - INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))))250 * NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)))) + snicind(ji) & 251 * (1 - NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji))))) 252 252 zdeltah(ji) = zdeltah(ji) + zh_s(ji) 253 253 END DO … … 258 258 ! 0 if not 259 259 DO ji = kideb, kiut 260 snicswi(ji) = MAX(0, INT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji)))))260 snicswi(ji) = MAX(0,NINT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji))))) 261 261 ENDDO 262 262 … … 279 279 280 280 DO ji = kideb, kiut 281 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1 .- snicind(ji) ) * snicswi(ji)281 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1 - snicind(ji) ) * snicswi(ji) 282 282 ! cotes of the top of the layers 283 283 zm0(ji,0) = 0._wp … … 291 291 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 292 292 limsum = MIN( limsum , nlay_s ) 293 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * limsum294 END DO 295 END DO 296 297 DO ji = kideb, kiut 298 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - snicswi(ji) * dh_snowice(ji) + zh_s(ji) * nlays0299 zm0(ji,1) = dh_s_tot(ji) * (1 -snswi(ji) ) + snswi(ji) * zm0(ji,1)293 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * REAL( limsum ) 294 END DO 295 END DO 296 297 DO ji = kideb, kiut 298 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - REAL( snicswi(ji) ) * dh_snowice(ji) + zh_s(ji) * REAL( nlays0 ) 299 zm0(ji,1) = dh_s_tot(ji) * REAL( 1 - snswi(ji) ) + REAL( snswi(ji) ) * zm0(ji,1) 300 300 END DO 301 301 … … 309 309 310 310 DO ji = kideb, kiut ! layer heat content 311 qm0 (ji,1) = rhosn * ( cpic * ( rtt - ( 1.- snswi(ji) ) * tatm_ice_1d(ji) &312 & - snswi(ji)* t_s_b (ji,1) ) &311 qm0 (ji,1) = rhosn * ( cpic * ( rtt - REAL( 1 - snswi(ji) ) * tatm_ice_1d(ji) & 312 & - REAL( snswi(ji) ) * t_s_b (ji,1) ) & 313 313 & + lfus ) * zthick0(ji,1) 314 314 zqts_in(ji) = zqts_in(ji) + qm0(ji,1) … … 320 320 limsum = MIN( limsum , nlay_s ) 321 321 qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk) 322 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, epsi20- ht_s_b(ji) ) )323 zqts_in(ji) = zqts_in(ji) + ( 1.- snswi(ji) ) * qm0(ji,jk) * zswitch322 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, - ht_s_b(ji) ) ) 323 zqts_in(ji) = zqts_in(ji) + REAL( 1 - snswi(ji) ) * qm0(ji,jk) * zswitch 324 324 END DO ! jk 325 325 END DO ! ji … … 360 360 !------------------- 361 361 DO ji = kideb, kiut 362 zh_s(ji) = ht_s_b(ji) / nlay_s362 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 363 363 z_s(ji,0) = 0._wp 364 364 ENDDO … … 366 366 DO jk = 1, nlay_s 367 367 DO ji = kideb, kiut 368 z_s(ji,jk) = zh_s(ji) * jk368 z_s(ji,jk) = zh_s(ji) * REAL( jk ) 369 369 END DO 370 370 END DO … … 394 394 & - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10)) 395 395 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) & 396 & * MAX(0.0,SIGN(1.0, nbot0(ji)-layer0+epsi20))396 & * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 397 397 END DO 398 398 END DO … … 410 410 DO ji = kideb, kiut 411 411 IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN 412 zji = MOD( npb(ji) - 1, jpi ) + 1 413 zjj = ( npb(ji) - 1 ) / jpi + 1 414 WRITE(numout,*) ' violation of heat conservation : ', & 415 ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 416 WRITE(numout,*) ' ji, jj : ', zji, zjj 412 ii = MOD( npb(ji) - 1, jpi ) + 1 413 ij = ( npb(ji) - 1 ) / jpi + 1 414 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 415 WRITE(numout,*) ' ji, jj : ', ii, ij 417 416 WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) 418 417 WRITE(numout,*) ' zqts_in : ', zqts_in (ji) * r1_rdtice … … 441 440 DO jk = 1, nlay_s 442 441 DO ji = kideb, kiut 443 zswitch = MAX ( 0.0 , SIGN ( 1.0, epsi20- ht_s_b(ji) ) )442 zswitch = MAX ( 0.0 , SIGN ( 1.0, - ht_s_b(ji) ) ) 444 443 t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 ) 445 444 END DO … … 480 479 limsum = ( (icsuswi(ji)*(icsuind(ji)+jk-1) + & 481 480 (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) 482 zm0(ji,jk)= icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) &483 + limsum* zh_i(ji)484 END DO 485 END DO 486 487 DO ji = kideb, kiut 488 zm0(ji,nbot0(ji)) = icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) + dh_i_bott(ji) &489 + zh_i(ji) * nlayi0490 zm0(ji,1) = snicswi(ji)*dh_snowice(ji) +(1-snicswi(ji))*zm0(ji,1)481 zm0(ji,jk)= REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) & 482 + REAL(limsum) * zh_i(ji) 483 END DO 484 END DO 485 486 DO ji = kideb, kiut 487 zm0(ji,nbot0(ji)) = REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) + dh_i_bott(ji) & 488 + zh_i(ji) * REAL(nlayi0) 489 zm0(ji,1) = REAL(snicswi(ji))*dh_snowice(ji) + REAL(1-snicswi(ji))*zm0(ji,1) 491 490 END DO 492 491 … … 521 520 !---------------------------- 522 521 DO ji = kideb, kiut 523 ztmelts = ( 1.0- icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice524 & + icboswi(ji)* (-tmut * s_i_new(ji) ) & ! case of forming ice522 ztmelts = REAL( 1 - icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice 523 & + REAL( icboswi(ji) ) * (-tmut * s_i_new(ji) ) & ! case of forming ice 525 524 & + rtt ! in Kelvin 526 525 … … 528 527 ztform = t_i_b(ji,nlay_i) 529 528 IF( num_sal == 2 ) ztform = t_bo_b(ji) 530 qm0(ji,nbot0(ji)) = ( 1.0- icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice531 & + icboswi(ji) * rhoic * ( cpic*(ztmelts-ztform) & ! case of forming ice529 qm0(ji,nbot0(ji)) = REAL( 1 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice 530 & + REAL( icboswi(ji) ) * rhoic * ( cpic*(ztmelts-ztform) & ! case of forming ice 532 531 + lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) ) & 533 532 - rcp*(ztmelts-rtt) ) * zthick0(ji,nbot0(ji) ) … … 540 539 ! energy of the flooding seawater 541 540 zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * & 542 (rhoic - rhosn) / rhoic * snicswi(ji) ! generally positive541 (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive 543 542 ! Heat conservation diagnostic 544 543 qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic … … 549 548 ! = enthalpy of snow + enthalpy of frozen water 550 549 zqsnic = zqsnow(ji) + zqsnic 551 qm0(ji,1) = snicswi(ji) * zqsnic +( 1 - snicswi(ji) ) * qm0(ji,1)550 qm0(ji,1) = REAL(snicswi(ji)) * zqsnic + REAL( 1 - snicswi(ji) ) * qm0(ji,1) 552 551 553 552 END DO ! ji … … 556 555 DO ji = kideb, kiut 557 556 ! Heat conservation 558 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06 +epsi20) ) &559 & * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + epsi20) )557 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06) ) & 558 & * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) 560 559 END DO 561 560 END DO … … 575 574 !------------------ 576 575 DO ji = kideb, kiut 577 zh_i(ji) = ht_i_b(ji) / nlay_i576 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 578 577 ENDDO 579 578 … … 606 605 q_i_b(ji,layer1) = q_i_b(ji,layer1) & 607 606 + zrl01(layer1,layer0)*qm0(ji,layer0) & 608 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06 +epsi20)) &609 * MAX(0.0,SIGN(1.0, nbot0(ji)-layer0+epsi20))607 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06)) & 608 * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 610 609 END DO 611 610 END DO … … 622 621 END DO 623 622 ! 624 DO ji = kideb, kiut 625 IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN 626 zji = MOD( npb(ji) - 1, jpi ) + 1 627 zjj = ( npb(ji) - 1 ) / jpi + 1 628 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 629 WRITE(numout,*) ' ji, jj : ', zji, zjj 630 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 631 WRITE(numout,*) ' zqti_in : ', zqti_in (ji) * r1_rdtice 632 WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice 633 WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) 634 WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) 635 WRITE(numout,*) ' dh_snowice:', dh_snowice(ji) 636 WRITE(numout,*) ' icsuswi : ', icsuswi(ji) 637 WRITE(numout,*) ' icboswi : ', icboswi(ji) 638 WRITE(numout,*) ' snicswi : ', snicswi(ji) 639 ENDIF 640 END DO 623 IF ( con_i ) THEN 624 DO ji = kideb, kiut 625 IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN 626 ii = MOD( npb(ji) - 1, jpi ) + 1 627 ij = ( npb(ji) - 1 ) / jpi + 1 628 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 629 WRITE(numout,*) ' ji, jj : ', ii, ij 630 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 631 WRITE(numout,*) ' zqti_in : ', zqti_in (ji) * r1_rdtice 632 WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice 633 WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) 634 WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) 635 WRITE(numout,*) ' dh_snowice:', dh_snowice(ji) 636 WRITE(numout,*) ' icsuswi : ', icsuswi(ji) 637 WRITE(numout,*) ' icboswi : ', icboswi(ji) 638 WRITE(numout,*) ' snicswi : ', snicswi(ji) 639 ENDIF 640 END DO 641 ENDIF 641 642 642 643 !---------------------- -
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r3625 r4161 46 46 47 47 !!---------------------------------------------------------------------- 48 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)48 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 49 49 !! $Id$ 50 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 78 78 !! update ht_s_b, ht_i_b and tbif_1d(:,:) 79 79 !!------------------------------------------------------------------------ 80 INTEGER 81 INTEGER 82 INTEGER :: zji, zjj, iter ! - -83 REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde! local scalars80 INTEGER :: ji,jj,jk,jl,jm ! dummy loop indices 81 INTEGER :: layer, nbpac ! local integers 82 INTEGER :: ii, ij, iter ! - - 83 REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde ! local scalars 84 84 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 85 85 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 86 REAL(wp) :: zcoef ! - -87 86 LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness 88 87 CHARACTER (len = 15) :: fieldid … … 160 159 DO ji = 1, jpi 161 160 !Energy of melting q(S,T) [J.m-3] 162 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * nlay_i161 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * REAL( nlay_i ) 163 162 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) ) ) !0 if no ice and 1 if yes 164 163 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb … … 342 341 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 343 342 DO ji = 1, nbpac 344 zji = MOD( npac(ji) - 1 , jpi ) + 1345 zjj = ( npac(ji) - 1 ) / jpi + 1346 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m( zji,zjj) )343 ii = MOD( npac(ji) - 1 , jpi ) + 1 344 ij = ( npac(ji) - 1 ) / jpi + 1 345 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij) ) 347 346 END DO 348 347 CASE ( 3 ) ! Sice = F(z) [multiyear ice] … … 389 388 END DO 390 389 391 !---------------------------------392 ! Salt flux due to new ice growth393 !---------------------------------394 ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine)395 DO ji = 1, nbpac396 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice397 rdm_ice_1d(ji) = rdm_ice_1d(ji) + rhoic * zv_newice(ji)398 END DO ! ji399 400 390 !------------------------------------ 401 391 ! Diags for energy conservation test 402 392 !------------------------------------ 403 393 DO ji = 1, nbpac 404 zji = MOD( npac(ji) - 1 , jpi ) + 1405 zjj = ( npac(ji) - 1 ) / jpi + 1394 ii = MOD( npac(ji) - 1 , jpi ) + 1 395 ij = ( npac(ji) - 1 ) / jpi + 1 406 396 ! 407 zde = ze_newice(ji) / unit_fac * area( zji,zjj) * zv_newice(ji)397 zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji) 408 398 ! 409