- Timestamp:
- 2015-07-21T13:25:36+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5146 r5621 24 24 USE wrk_nemo ! work arrays 25 25 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 26 USE sbc_oce, ONLY : lk_cpl27 26 28 27 IMPLICIT NONE … … 170 169 CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 171 170 CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe ) 172 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)173 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0)174 CALL wrk_alloc( jpij, 175 CALL wrk_alloc( jpij, nlay_i+3,3, ztrid )171 CALL wrk_alloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0 ) 172 CALL wrk_alloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0 ) 173 CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis ) 174 CALL wrk_alloc( jpij,nlay_i+3,3, ztrid ) 176 175 177 176 CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) … … 283 282 END DO 284 283 285 !286 284 !------------------------------------------------------------------------------| 287 285 ! 3) Iterative procedure begins | … … 679 677 END DO 680 678 681 DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1679 DO numeq = nlay_i + nlay_s, nlay_s + 2, -1 682 680 DO ji = kideb , kiut 683 681 jk = numeq - nlay_s - 1 … … 746 744 !-------------------------------------------------------------------------! 747 745 DO ji = kideb, kiut 748 ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)749 IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) )750 746 ! ! surface ice conduction flux 751 747 isnow(ji) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) … … 760 756 761 757 ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 762 IF ( ln_it_qnsice ) hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - ( qns_ice_1d(:) - zqns_ice_b(:) ) * a_i_1d(:) 758 IF ( ln_it_qnsice ) THEN 759 DO ji = kideb, kiut 760 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 761 END DO 762 END IF 763 763 764 764 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! … … 772 772 ENDIF 773 773 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 774 775 ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 776 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 774 777 END DO 775 776 hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - zhfx_err(:) * a_i_1d(:)777 778 778 779 !----------------------------------------- … … 787 788 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 788 789 ENDIF 789 END DO 790 791 ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m-2) 792 DO ji = kideb, kiut 793 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 794 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 795 END DO 796 790 ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 791 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji) 792 END DO 797 793 ! 798 794 CALL wrk_dealloc( jpij, numeqmin, numeqmax ) 799 795 CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 800 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 801 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 802 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 803 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 804 CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 805 CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 796 CALL wrk_dealloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe ) 797 CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 798 CALL wrk_dealloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 799 CALL wrk_dealloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis ) 800 CALL wrk_dealloc( jpij,nlay_i+3,3, ztrid ) 806 801 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 807 802 … … 824 819 DO jk = 1, nlay_i ! Sea ice energy of melting 825 820 DO ji = kideb, kiut 826 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 827 rswitch = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rt0) - epsi20 ) ) 828 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 829 & + lfus * ( 1.0 - rswitch * ( ztmelts-rt0 ) / MIN( t_i_1d(ji,jk) - rt0, -epsi20 ) ) & 830 & - rcp * ( ztmelts-rt0 ) ) 821 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 822 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point 823 ! (sometimes dif scheme produces abnormally high temperatures) 824 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 825 & + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) ) & 826 & - rcp * ( ztmelts-rt0 ) ) 831 827 END DO 832 828 END DO
Note: See TracChangeset
for help on using the changeset viewer.