- Timestamp:
- 2015-07-10T13:28:53+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
r4624 r5581 9 9 10 10 !!---------------------------------------------------------------------- 11 !! dia_hsb : Diagnose the conservation of ocean heat and salt contents, and volume 12 !! dia_hsb_rst : Read or write DIA file in restart file 13 !! dia_hsb_init : Initialization of the conservation diagnostic 14 !!---------------------------------------------------------------------- 11 15 USE oce ! ocean dynamics and tracers 12 16 USE dom_oce ! ocean space and time domain 13 17 USE phycst ! physical constants 14 18 USE sbc_oce ! surface thermohaline fluxes 15 USE in_out_manager ! I/O manager 19 USE sbcrnf ! river runoff 20 USE sbcisf ! ice shelves 16 21 USE domvvl ! vertical scale factors 17 22 USE traqsr ! penetrative solar radiation 18 23 USE trabbc ! bottom boundary condition 19 USE lib_mpp ! distributed memory computing library20 24 USE trabbc ! bottom boundary condition 21 25 USE bdy_par ! (for lk_bdy) 26 USE restart ! ocean restart 27 ! 28 USE iom ! I/O manager 29 USE in_out_manager ! I/O manager 30 USE lib_fortran ! glob_sum 31 USE lib_mpp ! distributed memory computing library 22 32 USE timing ! preformance summary 23 USE iom ! I/O manager 24 USE lib_fortran ! glob_sum 25 USE restart ! ocean restart 26 USE wrk_nemo ! work arrays 27 USE sbcrnf ! river runoffd 33 USE wrk_nemo ! work arrays 28 34 29 35 IMPLICIT NONE … … 36 42 LOGICAL, PUBLIC :: ln_diahsb !: check the heat and salt budgets 37 43 38 REAL(dp) :: surf_tot ! 39 REAL(dp) :: frc_t , frc_s , frc_v ! global forcing trends 40 REAL(dp) :: frc_wn_t , frc_wn_s ! global forcing trends 41 REAL(dp), DIMENSION(:,:) , ALLOCATABLE :: surf , ssh_ini ! 42 REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: hc_loc_ini, sc_loc_ini, e3t_ini ! 43 REAL(dp), DIMENSION(:,:) , ALLOCATABLE :: ssh_hc_loc_ini, ssh_sc_loc_ini 44 REAL(wp) :: surf_tot ! ocean surface 45 REAL(wp) :: frc_t, frc_s, frc_v ! global forcing trends 46 REAL(wp) :: frc_wn_t, frc_wn_s ! global forcing trends 47 ! 48 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf , ssh_ini ! 49 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: ssh_hc_loc_ini, ssh_sc_loc_ini ! 50 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: hc_loc_ini, sc_loc_ini, e3t_ini ! 44 51 45 52 !! * Substitutions 46 53 # include "domzgr_substitute.h90" 47 54 # include "vectopt_loop_substitute.h90" 48 49 55 !!---------------------------------------------------------------------- 50 56 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 52 58 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 53 59 !!---------------------------------------------------------------------- 54 55 60 CONTAINS 56 61 … … 67 72 !!--------------------------------------------------------------------------- 68 73 INTEGER, INTENT(in) :: kt ! ocean time-step index 69 !! 70 INTEGER :: jk ! dummy loop indice 71 REAL(dp) :: zdiff_hc , zdiff_sc ! heat and salt content variations 72 REAL(dp) :: zdiff_hc1 , zdiff_sc1 ! - - - - - - - - 73 REAL(dp) :: zdiff_v1 , zdiff_v2 ! volume variation 74 REAL(dp) :: zerr_hc1 , zerr_sc1 ! heat and salt content misfit 75 REAL(dp) :: zvol_tot ! volume 76 REAL(dp) :: z_frc_trd_t , z_frc_trd_s ! - - 77 REAL(dp) :: z_frc_trd_v ! - - 78 REAL(dp) :: z_wn_trd_t , z_wn_trd_s ! - - 79 REAL(dp) :: z_ssh_hc , z_ssh_sc ! - - 74 ! 75 INTEGER :: ji, jj, jk ! dummy loop indice 76 REAL(wp) :: zdiff_hc , zdiff_sc ! heat and salt content variations 77 REAL(wp) :: zdiff_hc1 , zdiff_sc1 ! - - - - 78 REAL(wp) :: zdiff_v1 , zdiff_v2 ! volume variation 79 REAL(wp) :: zerr_hc1 , zerr_sc1 ! heat and salt content misfit 80 REAL(wp) :: zvol_tot ! volume 81 REAL(wp) :: z_frc_trd_t , z_frc_trd_s ! - - 82 REAL(wp) :: z_frc_trd_v ! - - 83 REAL(wp) :: z_wn_trd_t , z_wn_trd_s ! - - 84 REAL(wp) :: z_ssh_hc , z_ssh_sc ! - - 85 REAL(wp), DIMENSION(:,:), POINTER :: z2d0, z2d1 80 86 !!--------------------------------------------------------------------------- 81 87 IF( nn_timing == 1 ) CALL timing_start('dia_hsb') 82 88 CALL wrk_alloc( jpi,jpj, z2d0, z2d1 ) 89 ! 90 tsn(:,:,:,1) = tsn(:,:,:,1) * tmask(:,:,:) ; tsb(:,:,:,1) = tsb(:,:,:,1) * tmask(:,:,:) ; 91 tsn(:,:,:,2) = tsn(:,:,:,2) * tmask(:,:,:) ; tsb(:,:,:,2) = tsb(:,:,:,2) * tmask(:,:,:) ; 83 92 ! ------------------------- ! 84 93 ! 1 - Trends due to forcing ! 85 94 ! ------------------------- ! 86 z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) ) ! volume fluxes87 z_frc_trd_t = glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) ) ! heat fluxes88 z_frc_trd_s = glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) ) ! salt fluxes89 ! Add runoff heat & salt input95 z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) ) * surf(:,:) ) ! volume fluxes 96 z_frc_trd_t = glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) ) ! heat fluxes 97 z_frc_trd_s = glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) ) ! salt fluxes 98 ! Add runoff heat & salt input 90 99 IF( ln_rnf ) z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) ) 91 100 IF( ln_rnf_sal) z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) ) 101 ! Add ice shelf heat & salt input 102 IF( nn_isf .GE. 1 ) THEN 103 z_frc_trd_t = z_frc_trd_t & 104 & + glob_sum( ( risf_tsc(:,:,jp_tem) - rdivisf * fwfisf(:,:) * (-1.9) * r1_rau0 ) * surf(:,:) ) 105 z_frc_trd_s = z_frc_trd_s + (1.0_wp - rdivisf) * glob_sum( risf_tsc(:,:,jp_sal) * surf(:,:) ) 106 ENDIF 92 107 93 108 ! Add penetrative solar radiation … … 97 112 ! 98 113 IF( .NOT. lk_vvl ) THEN 99 z_wn_trd_t = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) ) 100 z_wn_trd_s = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) ) 114 IF ( ln_isfcav ) THEN 115 DO ji=1,jpi 116 DO jj=1,jpj 117 z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem) 118 z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal) 119 ENDDO 120 ENDDO 121 ELSE 122 z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) 123 z2d1(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) 124 END IF 125 z_wn_trd_t = - glob_sum( z2d0 ) 126 z_wn_trd_s = - glob_sum( z2d1 ) 101 127 ENDIF 102 128 … … 113 139 ! 2 - Content variations ! 114 140 ! ------------------------ ! 115 zdiff_v2 = 0. d0116 zdiff_hc = 0. d0117 zdiff_sc = 0. d0141 zdiff_v2 = 0._wp 142 zdiff_hc = 0._wp 143 zdiff_sc = 0._wp 118 144 119 145 ! volume variation (calculated with ssh) … … 122 148 ! heat & salt content variation (associated with ssh) 123 149 IF( .NOT. lk_vvl ) THEN 124 z_ssh_hc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) ) 125 z_ssh_sc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) ) 150 IF ( ln_isfcav ) THEN 151 DO ji = 1, jpi 152 DO jj = 1, jpj 153 z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) ) 154 z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) ) 155 END DO 156 END DO 157 ELSE 158 z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) 159 z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) 160 END IF 161 z_ssh_hc = glob_sum( z2d0 ) 162 z_ssh_sc = glob_sum( z2d1 ) 126 163 ENDIF 127 164 … … 153 190 ! 3 - Diagnostics writing ! 154 191 ! ----------------------- ! 155 zvol_tot = 0.d0 ! total ocean volume192 zvol_tot = 0._wp ! total ocean volume (calculated with scale factors) 156 193 DO jk = 1, jpkm1 157 194 zvol_tot = zvol_tot + glob_sum( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) ) 158 195 END DO 196 197 !!gm to be added ? 198 ! IF( .NOT. lk_vvl ) THEN ! fixed volume, add the ssh contribution 199 ! zvol_tot = zvol_tot + glob_sum( surf(:,:) * sshn(:,:) ) 200 ! ENDIF 201 !!gm end 202 159 203 160 204 IF( lk_vvl ) THEN … … 183 227 IF( lrst_oce ) CALL dia_hsb_rst( kt, 'WRITE' ) 184 228 229 CALL wrk_dealloc( jpi,jpj, z2d0, z2d1 ) 230 185 231 IF( nn_timing == 1 ) CALL timing_stop('dia_hsb') 186 !232 ! 187 233 END SUBROUTINE dia_hsb 188 234 189 190 SUBROUTINE dia_hsb_init191 !!---------------------------------------------------------------------------192 !! *** ROUTINE dia_hsb ***193 !!194 !! ** Purpose: Initialization for the heat salt volume budgets195 !!196 !! ** Method : Compute initial heat content, salt content and volume197 !!198 !! ** Action : - Compute initial heat content, salt content and volume199 !! - Initialize forcing trends200 !! - Compute coefficients for conversion201 !!---------------------------------------------------------------------------202 INTEGER :: jk ! dummy loop indice203 INTEGER :: ierror ! local integer204 !!205 NAMELIST/namhsb/ ln_diahsb206 !207 INTEGER :: ios208 !!----------------------------------------------------------------------209 210 IF(lwp) THEN211 WRITE(numout,*)212 WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'213 WRITE(numout,*) '~~~~~~~~ '214 ENDIF215 216 REWIND( numnam_ref ) ! Namelist namhsb in reference namelist217 READ ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)218 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp )219 220 REWIND( numnam_cfg ) ! Namelist namhsb in configuration namelist221 READ ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )222 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp )223 IF(lwm) WRITE ( numond, namhsb )224 225 !226 IF(lwp) THEN ! Control print227 WRITE(numout,*)228 WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'229 WRITE(numout,*) '~~~~~~~~~~~~'230 WRITE(numout,*) ' Namelist namhsb : set hsb parameters'231 WRITE(numout,*) ' Switch for hsb diagnostic (T) or not (F) ln_diahsb = ', ln_diahsb232 WRITE(numout,*)233 ENDIF234 235 IF( .NOT. ln_diahsb ) RETURN236 ! IF( .NOT. lk_mpp_rep ) &237 ! CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', &238 ! & ' whereas the global sum to be precise must be done in double precision ',&239 ! & ' please add key_mpp_rep')240 241 ! ------------------- !242 ! 1 - Allocate memory !243 ! ------------------- !244 ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), &245 & e3t_ini(jpi,jpj,jpk), surf(jpi,jpj), ssh_ini(jpi,jpj), STAT=ierror )246 IF( ierror > 0 ) THEN247 CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) ; RETURN248 ENDIF249 250 IF(.NOT. lk_vvl ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )251 IF( ierror > 0 ) THEN252 CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) ; RETURN253 ENDIF254 255 ! ----------------------------------------------- !256 ! 2 - Time independant variables and file opening !257 ! ----------------------------------------------- !258 IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"259 IF(lwp) WRITE(numout,*) '~~~~~~~'260 surf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:) ! masked surface grid cell area261 surf_tot = glob_sum( surf(:,:) ) ! total ocean surface area262 263 IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )264 !265 ! ---------------------------------- !266 ! 4 - initial conservation variables !267 ! ---------------------------------- !268 CALL dia_hsb_rst( nit000, 'READ' ) !* read or initialize all required files269 !270 END SUBROUTINE dia_hsb_init271 235 272 236 SUBROUTINE dia_hsb_rst( kt, cdrw ) … … 281 245 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 282 246 ! 283 INTEGER :: j k !284 INTEGER :: id1 ! local integers247 INTEGER :: ji, jj, jk ! dummy loop indices 248 INTEGER :: id1 ! local integers 285 249 !!---------------------------------------------------------------------- 286 250 ! … … 317 281 sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) ! initial salt content 318 282 END DO 319 frc_v = 0. d0! volume trend due to forcing320 frc_t = 0. d0! heat content - - - -321 frc_s = 0. d0! salt content - - - -283 frc_v = 0._wp ! volume trend due to forcing 284 frc_t = 0._wp ! heat content - - - - 285 frc_s = 0._wp ! salt content - - - - 322 286 IF( .NOT. lk_vvl ) THEN 323 ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) ! initial heat content in ssh 324 ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:) ! initial salt content in ssh 325 frc_wn_t = 0.d0 ! initial heat content misfit due to free surface 326 frc_wn_s = 0.d0 ! initial salt content misfit due to free surface 287 IF ( ln_isfcav ) THEN 288 DO ji=1,jpi 289 DO jj=1,jpj 290 ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) ! initial heat content in ssh 291 ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) ! initial salt content in ssh 292 ENDDO 293 ENDDO 294 ELSE 295 ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) ! initial heat content in ssh 296 ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:) ! initial salt content in ssh 297 END IF 298 frc_wn_t = 0._wp ! initial heat content misfit due to free surface 299 frc_wn_s = 0._wp ! initial salt content misfit due to free surface 327 300 ENDIF 328 301 ENDIF … … 354 327 END SUBROUTINE dia_hsb_rst 355 328 329 330 SUBROUTINE dia_hsb_init 331 !!--------------------------------------------------------------------------- 332 !! *** ROUTINE dia_hsb *** 333 !! 334 !! ** Purpose: Initialization for the heat salt volume budgets 335 !! 336 !! ** Method : Compute initial heat content, salt content and volume 337 !! 338 !! ** Action : - Compute initial heat content, salt content and volume 339 !! - Initialize forcing trends 340 !! - Compute coefficients for conversion 341 !!--------------------------------------------------------------------------- 342 INTEGER :: jk ! dummy loop indice 343 INTEGER :: ierror ! local integer 344 INTEGER :: ios 345 ! 346 NAMELIST/namhsb/ ln_diahsb 347 !!---------------------------------------------------------------------- 348 349 IF(lwp) THEN 350 WRITE(numout,*) 351 WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets' 352 WRITE(numout,*) '~~~~~~~~ ' 353 ENDIF 354 355 REWIND( numnam_ref ) ! Namelist namhsb in reference namelist 356 READ ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901) 357 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp ) 358 359 REWIND( numnam_cfg ) ! Namelist namhsb in configuration namelist 360 READ ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 ) 361 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp ) 362 IF(lwm) WRITE ( numond, namhsb ) 363 364 ! 365 IF(lwp) THEN ! Control print 366 WRITE(numout,*) 367 WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets' 368 WRITE(numout,*) '~~~~~~~~~~~~' 369 WRITE(numout,*) ' Namelist namhsb : set hsb parameters' 370 WRITE(numout,*) ' Switch for hsb diagnostic (T) or not (F) ln_diahsb = ', ln_diahsb 371 WRITE(numout,*) 372 ENDIF 373 374 IF( .NOT. ln_diahsb ) RETURN 375 ! IF( .NOT. lk_mpp_rep ) & 376 ! CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', & 377 ! & ' whereas the global sum to be precise must be done in double precision ',& 378 ! & ' please add key_mpp_rep') 379 380 ! ------------------- ! 381 ! 1 - Allocate memory ! 382 ! ------------------- ! 383 ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), & 384 & e3t_ini(jpi,jpj,jpk), surf(jpi,jpj), ssh_ini(jpi,jpj), STAT=ierror ) 385 IF( ierror > 0 ) THEN 386 CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) ; RETURN 387 ENDIF 388 389 IF(.NOT. lk_vvl ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror ) 390 IF( ierror > 0 ) THEN 391 CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) ; RETURN 392 ENDIF 393 394 ! ----------------------------------------------- ! 395 ! 2 - Time independant variables and file opening ! 396 ! ----------------------------------------------- ! 397 IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated" 398 IF(lwp) WRITE(numout,*) '~~~~~~~' 399 surf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:) ! masked surface grid cell area 400 surf_tot = glob_sum( surf(:,:) ) ! total ocean surface area 401 402 IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' ) 403 ! 404 ! ---------------------------------- ! 405 ! 4 - initial conservation variables ! 406 ! ---------------------------------- ! 407 CALL dia_hsb_rst( nit000, 'READ' ) !* read or initialize all required files 408 ! 409 END SUBROUTINE dia_hsb_init 410 356 411 !!====================================================================== 357 412 END MODULE diahsb
Note: See TracChangeset
for help on using the changeset viewer.