Changeset 5989 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
- Timestamp:
- 2015-12-03T09:10:32+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r5260 r5989 37 37 USE limdyn ! Ice dynamics 38 38 USE limtrp ! Ice transport 39 USE limhdf ! Ice horizontal diffusion 39 40 USE limthd ! Ice thermodynamics 40 41 USE limitd_me ! Mechanics on ice thickness distribution … … 110 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_os, zalb_cs ! ice albedo under overcast/clear sky 111 112 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_ice ! mean ice albedo (for coupled) 113 REAL(wp), POINTER, DIMENSION(:,: ) :: zutau_ice, zvtau_ice 112 114 !!---------------------------------------------------------------------- 113 115 114 116 IF( nn_timing == 1 ) CALL timing_start('sbc_ice_lim') 115 117 116 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! Ice time-step only 117 !-----------------------! 118 ! --- Bulk Formulae --- ! 119 !-----------------------! 120 u_oce(:,:) = ssu_m(:,:) * umask(:,:,1) ! mean surface ocean current at ice velocity point 121 v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1) ! (C-grid dynamics : U- & V-points as the ocean) 118 !-----------------------! 119 ! --- Ice time step --- ! 120 !-----------------------! 121 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 122 123 ! mean surface ocean current at ice velocity point (C-grid dynamics : U- & V-points as the ocean) 124 u_oce(:,:) = ssu_m(:,:) * umask(:,:,1) 125 v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1) 122 126 123 127 ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 124 t_bo(:,:) = ( eos_fzp( sss_m ) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 125 ! 126 ! Ice albedo 127 CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 128 CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 129 130 ! CORE and COUPLED bulk formulations 131 SELECT CASE( kblk ) 132 CASE( jp_core , jp_cpl ) 133 134 ! albedo depends on cloud fraction because of non-linear spectral effects 135 zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 136 ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 137 ! (zalb_ice) is computed within the bulk routine 138 139 END SELECT 140 128 CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) 129 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 130 141 131 ! Mask sea ice surface temperature (set to rt0 over land) 142 132 DO jl = 1, jpl 143 133 t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 144 END DO 145 146 ! Bulk formulae - provides the following fields: 147 ! utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2] 134 END DO 135 ! 136 !------------------------------------------------! 137 ! --- Dynamical coupling with the atmosphere --- ! 138 !------------------------------------------------! 139 ! It provides the following fields: 140 ! utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2] 141 !----------------------------------------------------------------- 142 SELECT CASE( kblk ) 143 CASE( jp_clio ) ; CALL blk_ice_clio_tau ! CLIO bulk formulation 144 CASE( jp_core ) ; CALL blk_ice_core_tau ! CORE bulk formulation 145 CASE( jp_purecpl ) ; CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) ! Coupled formulation 146 END SELECT 147 148 IF( ln_mixcpl) THEN ! Case of a mixed Bulk/Coupled formulation 149 CALL wrk_alloc( jpi,jpj , zutau_ice, zvtau_ice) 150 CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice ) 151 utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 152 vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 153 CALL wrk_dealloc( jpi,jpj , zutau_ice, zvtau_ice) 154 ENDIF 155 156 !-------------------------------------------------------! 157 ! --- ice dynamics and transport (except in 1D case) ---! 158 !-------------------------------------------------------! 159 numit = numit + nn_fsbc ! Ice model time step 160 ! 161 CALL sbc_lim_bef ! Store previous ice values 162 CALL sbc_lim_diag0 ! set diag of mass, heat and salt fluxes to 0 163 CALL lim_rst_opn( kt ) ! Open Ice restart file 164 ! 165 IF( .NOT. lk_c1d ) THEN 166 ! 167 CALL lim_dyn( kt ) ! Ice dynamics ( rheology/dynamics ) 168 ! 169 CALL lim_trp( kt ) ! Ice transport ( Advection/diffusion ) 170 ! 171 IF( nn_monocat /= 2 ) CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting) 172 ! 173 #if defined key_bdy 174 CALL bdy_ice_lim( kt ) ! bdy ice thermo 175 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 176 #endif 177 ! 178 CALL lim_update1( kt ) ! Corrections 179 ! 180 ENDIF 181 182 ! previous lead fraction and ice volume for flux calculations 183 CALL sbc_lim_bef 184 CALL lim_var_glo2eqv ! ht_i and ht_s for ice albedo calculation 185 CALL lim_var_agg(1) ! at_i for coupling (via pfrld) 186 pfrld(:,:) = 1._wp - at_i(:,:) 187 phicif(:,:) = vt_i(:,:) 188 189 !------------------------------------------------------! 190 ! --- Thermodynamical coupling with the atmosphere --- ! 191 !------------------------------------------------------! 192 ! It provides the following fields: 148 193 ! qsr_ice , qns_ice : solar & non solar heat flux over ice (T-point) [W/m2] 149 194 ! qla_ice : latent heat flux over ice (T-point) [W/m2] … … 151 196 ! tprecip , sprecip : total & solid precipitation (T-point) [Kg/m2/s] 152 197 ! fr1_i0 , fr2_i0 : 1sr & 2nd fraction of qsr penetration in ice [%] 153 ! 198 !---------------------------------------------------------------------------------------- 199 CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 200 CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 201 154 202 SELECT CASE( kblk ) 155 203 CASE( jp_clio ) ! CLIO bulk formulation 156 CALL blk_ice_clio( t_su , zalb_cs , zalb_os , zalb_ice , & 157 & utau_ice , vtau_ice , qns_ice , qsr_ice , & 158 & qla_ice , dqns_ice , dqla_ice , & 159 & tprecip , sprecip , & 160 & fr1_i0 , fr2_i0 , cp_ice_msh, jpl ) 161 ! 162 IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice , & 163 & dqns_ice, qla_ice, dqla_ice, nn_limflx ) 164 204 ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 205 ! (zalb_ice) is computed within the bulk routine 206 CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, zalb_ice ) 207 IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 208 IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 165 209 CASE( jp_core ) ! CORE bulk formulation 166 CALL blk_ice_core( t_su , u_ice , v_ice , zalb_ice , &167 & utau_ice , vtau_ice , qns_ice , qsr_ice , &168 & qla_ice , dqns_ice , dqla_ice , &169 & tprecip , sprecip , &170 & fr1_i0 , fr2_i0 , cp_ice_msh, jpl)171 !172 IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice , &173 & dqns_ice, qla_ice, dqla_ice, nn_limflx)174 !175 CASE ( jp_cpl )176 177 CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )178 210 ! albedo depends on cloud fraction because of non-linear spectral effects 211 zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 212 CALL blk_ice_core_flx( t_su, zalb_ice ) 213 IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 214 IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 215 CASE ( jp_purecpl ) 216 ! albedo depends on cloud fraction because of non-linear spectral effects 217 zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 218 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 219 ! clem: evap_ice is forced to 0 in coupled mode for now 220 ! but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models 221 evap_ice (:,:,:) = 0._wp ; devap_ice (:,:,:) = 0._wp 222 IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 179 223 END SELECT 180 181 !------------------------------! 182 ! --- LIM-3 main time-step --- ! 183 !------------------------------! 184 numit = numit + nn_fsbc ! Ice model time step 185 ! 186 CALL sbc_lim_bef ! Store previous ice values 187 188 CALL sbc_lim_diag0 ! set diag of mass, heat and salt fluxes to 0 189 190 CALL lim_rst_opn( kt ) ! Open Ice restart file 191 ! 192 ! ---------------------------------------------- 193 ! ice dynamics and transport (except in 1D case) 194 ! ---------------------------------------------- 195 IF( .NOT. lk_c1d ) THEN 196 197 CALL lim_dyn( kt ) ! Ice dynamics ( rheology/dynamics ) 198 199 CALL lim_trp( kt ) ! Ice transport ( Advection/diffusion ) 200 201 IF( nn_monocat /= 2 ) CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting) 202 203 #if defined key_bdy 204 CALL bdy_ice_lim( kt ) ! bdy ice thermo 205 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 206 #endif 207 CALL lim_update1( kt ) 208 209 ENDIF 210 211 CALL sbc_lim_bef ! Store previous ice values 212 213 ! ---------------------------------------------- 214 ! ice thermodynamics 215 ! ---------------------------------------------- 216 CALL lim_var_agg(1) 217 218 ! previous lead fraction and ice volume for flux calculations 219 pfrld(:,:) = 1._wp - at_i(:,:) 220 phicif(:,:) = vt_i(:,:) 221 222 SELECT CASE( kblk ) 223 CASE ( jp_cpl ) 224 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 225 IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice , & 226 & dqns_ice, qla_ice, dqla_ice, nn_limflx ) 227 ! Latent heat flux is forced to 0 in coupled: it is included in qns (non-solar heat flux) 228 qla_ice (:,:,:) = 0._wp 229 dqla_ice (:,:,:) = 0._wp 230 END SELECT 231 ! 232 CALL lim_thd( kt ) ! Ice thermodynamics 233 224 CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 225 226 !----------------------------! 227 ! --- ice thermodynamics --- ! 228 !----------------------------! 229 CALL lim_thd( kt ) ! Ice thermodynamics 230 ! 234 231 CALL lim_update2( kt ) ! Corrections 235 232 ! … … 237 234 ! 238 235 IF(ln_limdiaout) CALL lim_diahsb ! Diagnostics and outputs 239 236 ! 240 237 CALL lim_wri( 1 ) ! Ice outputs 241 238 ! 242 239 IF( kt == nit000 .AND. ln_rstart ) & 243 240 & CALL iom_close( numrir ) ! close input ice restart file … … 247 244 IF( ln_icectl ) CALL lim_ctl( kt ) ! alerts in case of model crash 248 245 ! 249 CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )250 !251 246 ENDIF ! End sea-ice time step only 252 247 253 !--------------------------------! 254 ! --- at all ocean time step --- ! 255 !--------------------------------! 256 ! Update surface ocean stresses (only in ice-dynamic case) 257 ! otherwise the atm.-ocean stresses are used everywhere 248 !-------------------------! 249 ! --- Ocean time step --- ! 250 !-------------------------! 251 ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere 258 252 IF( ln_limdyn ) CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) ) ! using before instantaneous surf. currents 259 253 !!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!! 260 254 ! 261 IF( nn_timing == 1 ) 255 IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_lim') 262 256 ! 263 257 END SUBROUTINE sbc_ice_lim … … 300 294 ! 301 295 CALL lim_itd_init ! ice thickness distribution initialization 296 ! 297 CALL lim_hdf_init ! set ice horizontal diffusion computation parameters 302 298 ! 303 299 CALL lim_thd_init ! set ice thermodynics parameters … … 346 342 !!------------------------------------------------------------------- 347 343 INTEGER :: ios ! Local integer output status for namelist read 348 NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_ out,&344 NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir, & 349 345 & ln_limdyn, rn_amax, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt 350 346 !!------------------------------------------------------------------- … … 475 471 476 472 477 SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, & 478 & pdqn_ice, pqla_ice, pdql_ice, k_limflx ) 473 SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx ) 479 474 !!--------------------------------------------------------------------- 480 475 !! *** ROUTINE ice_lim_flx *** … … 494 489 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqsr_ice ! net solar flux 495 490 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdqn_ice ! non solar flux sensitivity 496 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p qla_ice ! latent heat flux497 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pd ql_ice ! latent heat fluxsensitivity491 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pevap_ice ! sublimation 492 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdevap_ice ! sublimation sensitivity 498 493 ! 499 494 INTEGER :: jl ! dummy loop index … … 504 499 REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m ! Mean solar heat flux over all categories 505 500 REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m ! Mean non solar heat flux over all categories 506 REAL(wp), POINTER, DIMENSION(:,:) :: z_ qla_m ! Mean latent heat fluxover all categories501 REAL(wp), POINTER, DIMENSION(:,:) :: z_evap_m ! Mean sublimation over all categories 507 502 REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m ! Mean d(qns)/dT over all categories 508 REAL(wp), POINTER, DIMENSION(:,:) :: z_d ql_m ! Mean d(qla)/dT over all categories503 REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories 509 504 !!---------------------------------------------------------------------- 510 505 … … 514 509 SELECT CASE( k_limflx ) !== averaged on all ice categories ==! 515 510 CASE( 0 , 1 ) 516 CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_ qla_m, z_dqn_m, z_dql_m)517 ! 518 z_qns_m (:,:) = fice_ice_ave ( pqns_ice (:,:,:) )519 z_qsr_m (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )520 z_dqn_m (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )521 z_ qla_m(:,:) = fice_ice_ave ( pqla_ice (:,:,:) )522 z_d ql_m(:,:) = fice_ice_ave ( pdql_ice (:,:,:) )511 CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m) 512 ! 513 z_qns_m (:,:) = fice_ice_ave ( pqns_ice (:,:,:) ) 514 z_qsr_m (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) ) 515 z_dqn_m (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) ) 516 z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) ) 517 z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) ) 523 518 DO jl = 1, jpl 524 pdqn_ice (:,:,jl) = z_dqn_m(:,:)525 pd ql_ice(:,:,jl) = z_dql_m(:,:)519 pdqn_ice (:,:,jl) = z_dqn_m(:,:) 520 pdevap_ice(:,:,jl) = z_devap_m(:,:) 526 521 END DO 527 522 ! 528 523 DO jl = 1, jpl 529 pqns_ice (:,:,jl) = z_qns_m(:,:)530 pqsr_ice (:,:,jl) = z_qsr_m(:,:)531 p qla_ice(:,:,jl) = z_qla_m(:,:)524 pqns_ice (:,:,jl) = z_qns_m(:,:) 525 pqsr_ice (:,:,jl) = z_qsr_m(:,:) 526 pevap_ice(:,:,jl) = z_evap_m(:,:) 532 527 END DO 533 528 ! 534 CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_ qla_m, z_dqn_m, z_dql_m)529 CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m) 535 530 END SELECT 536 531 … … 542 537 ztem_m(:,:) = fice_ice_ave ( ptn_ice (:,:,:) ) 543 538 DO jl = 1, jpl 544 pqns_ice (:,:,jl) = pqns_ice(:,:,jl) + pdqn_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))545 p qla_ice(:,:,jl) = pqla_ice(:,:,jl) + pdql_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))546 pqsr_ice (:,:,jl) = pqsr_ice(:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )539 pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) ) 540 pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) ) 541 pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) ) 547 542 END DO 548 543 ! … … 593 588 wfx_spr(:,:) = 0._wp ; 594 589 595 hfx_in (:,:) = 0._wp ; hfx_out(:,:) = 0._wp596 590 hfx_thd(:,:) = 0._wp ; 597 591 hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp … … 610 604 611 605 END SUBROUTINE sbc_lim_diag0 612 606 607 613 608 FUNCTION fice_cell_ave ( ptab ) 614 609 !!-------------------------------------------------------------------------- … … 620 615 621 616 fice_cell_ave (:,:) = 0.0_wp 622 623 617 DO jl = 1, jpl 624 618 fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
Note: See TracChangeset
for help on using the changeset viewer.
