Changeset 5630 for branches/UKMO/dev_r5107_restart_func_and_date/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
- Timestamp:
- 2015-07-23T18:05:51+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5107_restart_func_and_date/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r5500 r5630 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 128 t_bo(:,:) = ( eos_fzp( sss_m ) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 125 !126 ! Ice albedo127 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 albedos129 130 ! CORE and COUPLED bulk formulations131 SELECT CASE( kblk )132 CASE( jp_core , jp_cpl )133 134 ! albedo depends on cloud fraction because of non-linear spectral effects135 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 albedo137 ! (zalb_ice) is computed within the bulk routine138 139 END SELECT140 129 141 130 ! Mask sea ice surface temperature (set to rt0 over land) 142 131 DO jl = 1, jpl 143 132 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] 133 END DO 134 ! 135 !------------------------------------------------! 136 ! --- Dynamical coupling with the atmosphere --- ! 137 !------------------------------------------------! 138 ! It provides the following fields: 139 ! utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2] 140 !----------------------------------------------------------------- 141 SELECT CASE( kblk ) 142 CASE( jp_clio ) ; CALL blk_ice_clio_tau ! CLIO bulk formulation 143 CASE( jp_core ) ; CALL blk_ice_core_tau ! CORE bulk formulation 144 CASE( jp_purecpl ) ; CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) ! Coupled formulation 145 END SELECT 146 147 IF( ln_mixcpl) THEN ! Case of a mixed Bulk/Coupled formulation 148 CALL wrk_alloc( jpi,jpj , zutau_ice, zvtau_ice) 149 CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice ) 150 utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 151 vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 152 CALL wrk_dealloc( jpi,jpj , zutau_ice, zvtau_ice) 153 ENDIF 154 155 !-------------------------------------------------------! 156 ! --- ice dynamics and transport (except in 1D case) ---! 157 !-------------------------------------------------------! 158 numit = numit + nn_fsbc ! Ice model time step 159 ! 160 CALL sbc_lim_bef ! Store previous ice values 161 CALL sbc_lim_diag0 ! set diag of mass, heat and salt fluxes to 0 162 CALL lim_rst_opn( kt ) ! Open Ice restart file 163 ! 164 IF( .NOT. lk_c1d ) THEN 165 ! 166 CALL lim_dyn( kt ) ! Ice dynamics ( rheology/dynamics ) 167 ! 168 CALL lim_trp( kt ) ! Ice transport ( Advection/diffusion ) 169 ! 170 IF( nn_monocat /= 2 ) CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting) 171 ! 172 #if defined key_bdy 173 CALL bdy_ice_lim( kt ) ! bdy ice thermo 174 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 175 #endif 176 ! 177 CALL lim_update1( kt ) ! Corrections 178 ! 179 ENDIF 180 181 ! previous lead fraction and ice volume for flux calculations 182 CALL sbc_lim_bef 183 CALL lim_var_glo2eqv ! ht_i and ht_s for ice albedo calculation 184 CALL lim_var_agg(1) ! at_i for coupling (via pfrld) 185 pfrld(:,:) = 1._wp - at_i(:,:) 186 phicif(:,:) = vt_i(:,:) 187 188 !------------------------------------------------------! 189 ! --- Thermodynamical coupling with the atmosphere --- ! 190 !------------------------------------------------------! 191 ! It provides the following fields: 148 192 ! qsr_ice , qns_ice : solar & non solar heat flux over ice (T-point) [W/m2] 149 193 ! qla_ice : latent heat flux over ice (T-point) [W/m2] … … 151 195 ! tprecip , sprecip : total & solid precipitation (T-point) [Kg/m2/s] 152 196 ! fr1_i0 , fr2_i0 : 1sr & 2nd fraction of qsr penetration in ice [%] 153 ! 197 !---------------------------------------------------------------------------------------- 198 CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 199 CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 200 154 201 SELECT CASE( kblk ) 155 202 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 203 ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 204 ! (zalb_ice) is computed within the bulk routine 205 CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, zalb_ice ) 206 IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 207 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 208 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 209 ! albedo depends on cloud fraction because of non-linear spectral effects 210 zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 211 CALL blk_ice_core_flx( t_su, zalb_ice ) 212 IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 213 IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 214 CASE ( jp_purecpl ) 215 ! albedo depends on cloud fraction because of non-linear spectral effects 216 zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 217 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 218 ! clem: evap_ice is forced to 0 in coupled mode for now 219 ! but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models 220 evap_ice (:,:,:) = 0._wp ; devap_ice (:,:,:) = 0._wp 221 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 222 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 223 CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 224 225 !----------------------------! 226 ! --- ice thermodynamics --- ! 227 !----------------------------! 228 CALL lim_thd( kt ) ! Ice thermodynamics 229 ! 234 230 CALL lim_update2( kt ) ! Corrections 235 231 ! … … 237 233 ! 238 234 IF(ln_limdiaout) CALL lim_diahsb ! Diagnostics and outputs 239 235 ! 240 236 CALL lim_wri( 1 ) ! Ice outputs 241 237 ! 242 238 IF( kt == nit000 .AND. ln_rstart ) & 243 239 & CALL iom_close( numrir ) ! close input ice restart file … … 247 243 IF( ln_icectl ) CALL lim_ctl( kt ) ! alerts in case of model crash 248 244 ! 249 CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )250 !251 245 ENDIF ! End sea-ice time step only 252 246 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 247 !-------------------------! 248 ! --- Ocean time step --- ! 249 !-------------------------! 250 ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere 258 251 IF( ln_limdyn ) CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) ) ! using before instantaneous surf. currents 259 252 !!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!! 260 253 ! 261 IF( nn_timing == 1 ) 254 IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_lim') 262 255 ! 263 256 END SUBROUTINE sbc_ice_lim … … 300 293 ! 301 294 CALL lim_itd_init ! ice thickness distribution initialization 295 ! 296 CALL lim_hdf_init ! set ice horizontal diffusion computation parameters 302 297 ! 303 298 CALL lim_thd_init ! set ice thermodynics parameters … … 475 470 476 471 477 SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, & 478 & pdqn_ice, pqla_ice, pdql_ice, k_limflx ) 472 SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx ) 479 473 !!--------------------------------------------------------------------- 480 474 !! *** ROUTINE ice_lim_flx *** … … 494 488 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqsr_ice ! net solar flux 495 489 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 fluxsensitivity490 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pevap_ice ! sublimation 491 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdevap_ice ! sublimation sensitivity 498 492 ! 499 493 INTEGER :: jl ! dummy loop index … … 504 498 REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m ! Mean solar heat flux over all categories 505 499 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 categories500 REAL(wp), POINTER, DIMENSION(:,:) :: z_evap_m ! Mean sublimation over all categories 507 501 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 categories502 REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories 509 503 !!---------------------------------------------------------------------- 510 504 … … 514 508 SELECT CASE( k_limflx ) !== averaged on all ice categories ==! 515 509 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 (:,:,:) )510 CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m) 511 ! 512 z_qns_m (:,:) = fice_ice_ave ( pqns_ice (:,:,:) ) 513 z_qsr_m (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) ) 514 z_dqn_m (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) ) 515 z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) ) 516 z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) ) 523 517 DO jl = 1, jpl 524 pdqn_ice (:,:,jl) = z_dqn_m(:,:)525 pd ql_ice(:,:,jl) = z_dql_m(:,:)518 pdqn_ice (:,:,jl) = z_dqn_m(:,:) 519 pdevap_ice(:,:,jl) = z_devap_m(:,:) 526 520 END DO 527 521 ! 528 522 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(:,:)523 pqns_ice (:,:,jl) = z_qns_m(:,:) 524 pqsr_ice (:,:,jl) = z_qsr_m(:,:) 525 pevap_ice(:,:,jl) = z_evap_m(:,:) 532 526 END DO 533 527 ! 534 CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_ qla_m, z_dqn_m, z_dql_m)528 CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m) 535 529 END SELECT 536 530 … … 542 536 ztem_m(:,:) = fice_ice_ave ( ptn_ice (:,:,:) ) 543 537 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(:,:) )538 pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) ) 539 pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) ) 540 pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) ) 547 541 END DO 548 542 ! … … 593 587 wfx_spr(:,:) = 0._wp ; 594 588 595 hfx_in (:,:) = 0._wp ; hfx_out(:,:) = 0._wp596 589 hfx_thd(:,:) = 0._wp ; 597 590 hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp … … 610 603 611 604 END SUBROUTINE sbc_lim_diag0 612 605 606 613 607 FUNCTION fice_cell_ave ( ptab ) 614 608 !!-------------------------------------------------------------------------- … … 620 614 621 615 fice_cell_ave (:,:) = 0.0_wp 622 623 616 DO jl = 1, jpl 624 617 fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
Note: See TracChangeset
for help on using the changeset viewer.
