Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
r2471 r2528 4 4 !! Ocean fluxes : domain averaged freshwater budget 5 5 !!====================================================================== 6 !! History : 8.2 !01-02 (E. Durand) Original code7 !! 8.5 !02-06 (G. Madec) F90: Free form and module8 !! 9.0 !06-08 (G. Madec) Surface module9 !! 9.2 !09-07 (C. Talandier) emp mean s spread over erp area6 !! History : OPA ! 2001-02 (E. Durand) Original code 7 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 8 !! 3.0 ! 2006-08 (G. Madec) Surface module 9 !! 3.2 ! 2009-07 (C. Talandier) emp mean s spread over erp area 10 10 !!---------------------------------------------------------------------- 11 11 … … 23 23 USE lib_mpp ! distribued memory computing library 24 24 USE lbclnk ! ocean lateral boundary conditions 25 USE lib_fortran 25 26 26 27 IMPLICIT NONE … … 31 32 REAL(wp) :: a_fwb_b ! annual domain averaged freshwater budget 32 33 REAL(wp) :: a_fwb ! for 2 year before (_b) and before year. 33 REAL(wp) :: empold ! empold to be suppressed34 REAL(wp) :: fwfold ! fwfold to be suppressed 34 35 REAL(wp) :: area ! global mean ocean surface (interior domain) 35 36 36 REAL(wp), DIMENSION(jpi,jpj) :: e1e2 _i ! area of the interior domain (e1t*e2t*tmask_i)37 REAL(wp), DIMENSION(jpi,jpj) :: e1e2 ! area of the interior domain (e1t*e2t) 37 38 38 39 !! * Substitutions … … 40 41 # include "vectopt_loop_substitute.h90" 41 42 !!---------------------------------------------------------------------- 42 !! OPA 9.0 , LOCEAN-IPSL (2006)43 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 43 44 !! $Id$ 44 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)45 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 45 46 !!---------------------------------------------------------------------- 46 47 CONTAINS … … 65 66 INTEGER :: inum ! temporary logical unit 66 67 INTEGER :: ikty, iyear ! 67 REAL(wp) :: z_ emp, z_emp_nsrf, zsum_emp, zsum_erp ! temporary scalars68 REAL(wp) :: z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp ! temporary scalars 68 69 REAL(wp) :: zsurf_neg, zsurf_pos, zsurf_tospread 69 70 REAL(wp), DIMENSION(jpi,jpj) :: ztmsk_neg, ztmsk_pos, ztmsk_tospread … … 72 73 ! 73 74 IF( kt == nit000 ) THEN 74 !75 75 IF(lwp) THEN 76 76 WRITE(numout,*) … … 79 79 IF( kn_fwb == 1 ) WRITE(numout,*) ' instantaneously set to zero' 80 80 IF( kn_fwb == 2 ) WRITE(numout,*) ' adjusted from previous year budget' 81 IF( kn_fwb == 3 ) WRITE(numout,*) ' empset to zero and spread out over erp area'81 IF( kn_fwb == 3 ) WRITE(numout,*) ' fwf set to zero and spread out over erp area' 82 82 ENDIF 83 83 ! 84 84 IF( kn_fwb == 3 .AND. nn_sssr /= 2 ) CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' ) 85 85 ! 86 e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:) 87 area = SUM( e1e2_i(:,:) ) 88 IF( lk_mpp ) CALL mpp_sum( area ) ! sum over the global domain 89 ! 86 e1e2(:,:) = e1t(:,:) * e2t(:,:) 87 area = glob_sum( e1e2(:,:) ) ! interior global domain surface 90 88 ENDIF 91 89 … … 93 91 SELECT CASE ( kn_fwb ) 94 92 ! 95 CASE ( 0 ) 96 WRITE(ctmp1,*)'sbc_fwb : nn_fwb = ', kn_fwb, ' is not yet associated to an option, choose either 1/2' 97 CALL ctl_stop( ctmp1 ) 93 CASE ( 1 ) !== global mean fwf set to zero ==! 98 94 ! 99 100 !101 CASE ( 1 ) ! global mean emp set to zero102 95 IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 103 z_emp = SUM( e1e2_i(:,:) * emp(:,:) ) / area 104 IF( lk_mpp ) CALL mpp_sum( z_emp ) ! sum over the global domain 105 emp (:,:) = emp (:,:) - z_emp 106 emps(:,:) = emps(:,:) - z_emp 96 z_fwf = glob_sum( e1e2(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area ! sum over the global domain 97 emp (:,:) = emp (:,:) - z_fwf 98 emps(:,:) = emps(:,:) - z_fwf 107 99 ENDIF 108 100 ! 109 CASE ( 2 ) ! emp budget adjusted from the previous year110 ! initialisation111 IF( kt == nit000 ) THEN 112 ! Read the corrective factor on precipitations (empold)101 CASE ( 2 ) !== fwf budget adjusted from the previous year ==! 102 ! 103 IF( kt == nit000 ) THEN ! initialisation 104 ! ! Read the corrective factor on precipitations (fwfold) 113 105 CALL ctl_opn( inum, 'EMPave_old.dat', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 114 106 READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb 115 107 CLOSE( inum ) 116 empold = a_fwb! current year freshwater budget correction117 ! ! estimate from the previous year budget108 fwfold = a_fwb ! current year freshwater budget correction 109 ! ! estimate from the previous year budget 118 110 IF(lwp)WRITE(numout,*) 119 IF(lwp)WRITE(numout,*)'sbc_fwb : year = ',iyear , ' freshwater budget correction = ', empold111 IF(lwp)WRITE(numout,*)'sbc_fwb : year = ',iyear , ' freshwater budget correction = ', fwfold 120 112 IF(lwp)WRITE(numout,*)' year = ',iyear-1, ' freshwater budget read = ', a_fwb 121 113 IF(lwp)WRITE(numout,*)' year = ',iyear-2, ' freshwater budget read = ', a_fwb_b 122 114 ENDIF 123 ! 124 ! Update empold if new year start 115 ! ! Update fwfold if new year start 125 116 ikty = 365 * 86400 / rdttra(1) !!bug use of 365 days leap year or 360d year !!!!!!! 126 117 IF( MOD( kt, ikty ) == 0 ) THEN 127 118 a_fwb_b = a_fwb 128 a_fwb = SUM( e1e2_i(:,:) * sshn(:,:) ) 129 IF( lk_mpp ) CALL mpp_sum( a_fwb ) ! sum over the global domain 119 a_fwb = glob_sum( e1e2(:,:) * sshn(:,:) ) ! sum over the global domain 130 120 a_fwb = a_fwb * 1.e+3 / ( area * 86400. * 365. ) ! convert in Kg/m3/s = mm/s 131 121 !!gm ! !!bug 365d year 132 empold = a_fwb! current year freshwater budget correction133 ! ! estimate from the previous year budget122 fwfold = a_fwb ! current year freshwater budget correction 123 ! ! estimate from the previous year budget 134 124 ENDIF 135 125 ! 136 ! correct the freshwater fluxes 137 IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 138 emp (:,:) = emp (:,:) + empold 139 emps(:,:) = emps(:,:) + empold 126 IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! correct the freshwater fluxes 127 emp (:,:) = emp (:,:) + fwfold 128 emps(:,:) = emps(:,:) + fwfold 140 129 ENDIF 141 130 ! 142 ! save empold value in a file 143 IF( kt == nitend .AND. lwp ) THEN 131 IF( kt == nitend .AND. lwp ) THEN ! save fwfold value in a file 144 132 CALL ctl_opn( inum, 'EMPave.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea ) 145 133 WRITE( inum, "(24X,I8,2ES24.16)" ) nyear, a_fwb_b, a_fwb … … 147 135 ENDIF 148 136 ! 149 CASE ( 3 ) ! global emp set to zero and spread out over erp area137 CASE ( 3 ) !== global fwf set to zero and spread out over erp area ==! 150 138 ! 151 139 IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 152 ! Select <0 and >0 area of erp 153 ztmsk_pos(:,:) = tmask_i(:,:) 154 WHERE( erp < 0.e0 ) ztmsk_pos = 0.e0 140 ztmsk_pos(:,:) = tmask_i(:,:) ! Select <0 and >0 area of erp 141 WHERE( erp < 0._wp ) ztmsk_pos = 0._wp 155 142 ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:) 156 157 ! Area filled by <0 and >0 erp158 zsurf_neg = SUM( e1e2_i(:,:)*ztmsk_neg(:,:) )159 zsurf_pos = SUM( e1e2_i(:,:)*ztmsk_pos(:,:) )160 161 ! emp global mean162 z_emp = SUM( e1e2_i(:,:) * emp(:,:) ) / area163 143 ! 164 IF( lk_mpp ) CALL mpp_sum( z_emp )165 IF( lk_mpp ) CALL mpp_sum( zsurf_neg)166 IF( lk_mpp ) CALL mpp_sum( zsurf_pos )167 168 IF( z_emp < 0.e0 ) THEN169 ! to spread out over >0 erp area to increase evaporation damping process170 zsurf_tospread = zsurf_pos144 zsurf_neg = glob_sum( e1e2(:,:)*ztmsk_neg(:,:) ) ! Area filled by <0 and >0 erp 145 zsurf_pos = glob_sum( e1e2(:,:)*ztmsk_pos(:,:) ) 146 ! ! fwf global mean 147 z_fwf = glob_sum( e1e2(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area 148 ! 149 IF( z_fwf < 0._wp ) THEN ! spread out over >0 erp area to increase evaporation 150 zsurf_tospread = zsurf_pos 171 151 ztmsk_tospread(:,:) = ztmsk_pos(:,:) 172 ELSE 173 ! to spread out over <0 erp area to increase precipitation damping process 174 zsurf_tospread = zsurf_neg 152 ELSE ! spread out over <0 erp area to increase precipitation 153 zsurf_tospread = zsurf_neg 175 154 ztmsk_tospread(:,:) = ztmsk_neg(:,:) 176 155 ENDIF 177 178 ! emp global mean over <0 or >0 erp area 179 zsum_emp = SUM( e1e2_i(:,:) * z_emp ) 180 IF( lk_mpp ) CALL mpp_sum( zsum_emp ) 181 z_emp_nsrf = zsum_emp / ( zsurf_tospread + rsmall ) 182 ! weight to respect erp field 2D structure 183 zsum_erp = SUM( ztmsk_tospread(:,:) * erp(:,:) * e1e2_i(:,:) ) 184 IF( lk_mpp ) CALL mpp_sum( zsum_erp ) 156 ! 157 zsum_fwf = glob_sum( e1e2(:,:) * z_fwf ) ! fwf global mean over <0 or >0 erp area 158 !!gm : zsum_fwf = z_fwf * area ??? it is right? I think so.... 159 z_fwf_nsrf = zsum_fwf / ( zsurf_tospread + rsmall ) 160 ! ! weight to respect erp field 2D structure 161 zsum_erp = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2(:,:) ) 185 162 z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall ) 186 187 ! final correction term to apply188 zerp_cor(:,:) = -1. * z_emp_nsrf * zsurf_tospread * z_wgt(:,:)189 163 ! ! final correction term to apply 164 zerp_cor(:,:) = -1. * z_fwf_nsrf * zsurf_tospread * z_wgt(:,:) 165 ! 166 !!gm ===>>>> lbc_lnk should be useless as all the computation is done over the whole domain ! 190 167 CALL lbc_lnk( zerp_cor, 'T', 1. ) 191 168 ! 192 169 emp (:,:) = emp (:,:) + zerp_cor(:,:) 193 170 emps(:,:) = emps(:,:) + zerp_cor(:,:) 194 171 erp (:,:) = erp (:,:) + zerp_cor(:,:) 195 196 IF( nprint == 1 .AND. lwp ) THEN 197 IF( z_ emp < 0.e0) THEN198 WRITE(numout,*)' z_emp< 0'199 WRITE(numout,*)' SUM(erp+) = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2_i(:,:) )*1.e-3,' m3.s-1'172 ! 173 IF( nprint == 1 .AND. lwp ) THEN ! control print 174 IF( z_fwf < 0._wp ) THEN 175 WRITE(numout,*)' z_fwf < 0' 176 WRITE(numout,*)' SUM(erp+) = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2(:,:) )*1.e-9,' Sv' 200 177 ELSE 201 WRITE(numout,*)' z_emp>= 0'202 WRITE(numout,*)' SUM(erp-) = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2_i(:,:) )*1.e-3,' m3.s-1'178 WRITE(numout,*)' z_fwf >= 0' 179 WRITE(numout,*)' SUM(erp-) = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2(:,:) )*1.e-9,' Sv' 203 180 ENDIF 204 WRITE(numout,*)' SUM(empG) = ', SUM( z_emp*e1e2_i(:,:) )*1.e-3,' m3.s-1'205 WRITE(numout,*)' z_emp = ', z_emp ,' mm.s-1'206 WRITE(numout,*)' z_emp_nsrf = ', z_emp_nsrf ,' mm.s-1'207 WRITE(numout,*)' MIN(zerp_cor)= ', MINVAL(zerp_cor)208 WRITE(numout,*)' MAX(zerp_cor)= ', MAXVAL(zerp_cor)181 WRITE(numout,*)' SUM(empG) = ', SUM( z_fwf*e1e2(:,:) )*1.e-9,' Sv' 182 WRITE(numout,*)' z_fwf = ', z_fwf ,' Kg/m2/s' 183 WRITE(numout,*)' z_fwf_nsrf = ', z_fwf_nsrf ,' Kg/m2/s' 184 WRITE(numout,*)' MIN(zerp_cor) = ', MINVAL(zerp_cor) 185 WRITE(numout,*)' MAX(zerp_cor) = ', MAXVAL(zerp_cor) 209 186 ENDIF 210 !211 187 ENDIF 212 188 ! 213 CASE DEFAULT ! you should never be there 214 WRITE(ctmp1,*)'sbc_fwb : nn_fwb = ', kn_fwb, ' is not permitted for the FreshWater Budget correction, choose either 1/2' 215 CALL ctl_stop( ctmp1 ) 189 CASE DEFAULT !== you should never be there ==! 190 CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1, 2 or 3' ) 216 191 ! 217 192 END SELECT
Note: See TracChangeset
for help on using the changeset viewer.