Changeset 5350 for branches/2014/dev_r5134_UKMO4_CF_compliance/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
- Timestamp:
- 2015-06-04T16:12:19+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r5134_UKMO4_CF_compliance/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5128 r5350 120 120 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 121 121 REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface 122 REAL(wp), POINTER, DIMENSION(:) :: zqns_ice_b ! solar radiation absorbed at the surface 122 123 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 123 124 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function … … 168 169 CALL wrk_alloc( jpij, numeqmin, numeqmax ) 169 170 CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 170 CALL wrk_alloc( jpij, zf, dzf, z errit, zdifcase, zftrice, zihic, zghe )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, 174 CALL wrk_alloc( jpij, nlay_i+3,3, ztrid )171 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,nlay_i+3, zindterm, zindtbis, zdiagbis ) 175 CALL wrk_alloc( jpij,nlay_i+3,3, ztrid ) 175 176 176 177 CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) … … 242 243 !------------------------------------------------------- 243 244 DO ji = kideb , kiut 244 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 245 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 246 dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 245 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 246 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 247 dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 248 zqns_ice_b(ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value 247 249 END DO 248 250 … … 452 454 !------------------------------------------------------------------------------| 453 455 ! 454 IF ( .NOT. lk_cpl ) THEN !--- forced atmosphere case456 IF ( ln_it_qnsice ) THEN 455 457 DO ji = kideb , kiut 456 458 ! update of the non solar flux according to the update in T_su … … 677 679 END DO 678 680 679 DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1681 DO numeq = nlay_i + nlay_s, nlay_s + 2, -1 680 682 DO ji = kideb , kiut 681 683 jk = numeq - nlay_s - 1 … … 757 759 CALL lim_thd_enmelt( kideb, kiut ) 758 760 761 ! --- 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(:) 759 763 760 764 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! … … 768 772 ENDIF 769 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) 770 777 END DO 771 772 ! diagnose external surface (forced case) or bottom (forced case) from heat conservation773 IF( .NOT. lk_cpl ) THEN ! --- forced case: qns_ice and fc_su are diagnosed774 !775 DO ji = kideb, kiut776 qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err(ji)777 fc_su (ji) = fc_su(ji) - zhfx_err(ji)778 END DO779 !780 ELSE ! --- coupled case: ocean turbulent heat flux is diagnosed781 !782 DO ji = kideb, kiut783 fhtur_1d (ji) = fhtur_1d(ji) - zhfx_err(ji)784 END DO785 !786 ENDIF787 778 788 779 !----------------------------------------- … … 797 788 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 798 789 ENDIF 799 END DO 800 801 ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m-2) 802 DO ji = kideb, kiut 803 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 804 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 805 END DO 806 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 807 793 ! 808 794 CALL wrk_dealloc( jpij, numeqmin, numeqmax ) 809 795 CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 810 796 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 811 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 812 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 813 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 814 CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 815 CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 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 ) 816 801 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 817 802 … … 834 819 DO jk = 1, nlay_i ! Sea ice energy of melting 835 820 DO ji = kideb, kiut 836 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 837 rswitch = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rt0) - epsi20 ) ) 838 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 839 & + lfus * ( 1.0 - rswitch * ( ztmelts-rt0 ) / MIN( t_i_1d(ji,jk) - rt0, -epsi20 ) ) & 840 & - 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 ) ) 841 827 END DO 842 828 END DO
Note: See TracChangeset
for help on using the changeset viewer.