Changeset 13489 for branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1_v2/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90
- Timestamp:
- 2020-09-17T19:13:57+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1_v2/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90
r10728 r13489 5 5 !!====================================================================== 6 6 !! History : 1.0 ! 2015-04 (PML) Original code 7 !! History : 1.1 ! 2020-06 (PML) Update to FABM 1.0, improved performance 7 8 !!---------------------------------------------------------------------- 8 9 #if defined key_fabm … … 19 20 USE trd_oce 20 21 USE trdtrc 22 USE traqsr,only: ln_qsr_spec 21 23 #if defined key_trdtrc && defined key_iomput 22 24 USE trdtrc_oce … … 24 26 25 27 USE oce, only: tsn ! Needed? 26 USE sbc_oce, only: lk_oasis 28 USE sbc_oce, only: lk_oasis,fr_i 27 29 USE dom_oce 28 30 USE zdf_oce 29 !USE iom31 USE iom 30 32 USE xios 31 33 USE cpl_oasis3 … … 36 38 USE asmbgc, ONLY: mld_choice_bgc 37 39 USE lbclnk 40 USE fabm_types, ONLY: standard_variables 41 USE par_fabm, ONLY: jp_fabm_swr 38 42 39 43 !USE fldread ! time interpolation 40 41 USE fabm42 44 43 45 IMPLICIT NONE … … 48 50 PRIVATE 49 51 50 PUBLIC trc_sms_fabm ! called by trcsms.F90 module 51 PUBLIC trc_sms_fabm_alloc ! called by trcini_fabm.F90 module 52 PUBLIC trc_sms_fabm_check_mass 53 PUBLIC st2d_fabm_nxt ! 2D state intergration 54 PUBLIC compute_fabm ! Compute FABM sources, sinks and diagnostics 55 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: flux ! Cross-interface flux of pelagic variables (# m-2 s-1) 57 58 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ext ! Light extinction coefficient (m-1) 59 60 ! Work array for mass aggregation 61 REAL(wp), ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: current_total 62 63 64 ! Arrays for environmental variables 65 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: prn,rho 66 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: taubot 67 68 ! repair counters 69 INTEGER :: repair_interior_count,repair_surface_count,repair_bottom_count 70 71 ! state check type 72 TYPE type_state 73 LOGICAL :: valid 74 LOGICAL :: repaired 75 END TYPE 76 77 REAL(wp), PUBLIC :: daynumber_in_year 78 79 TYPE (type_bulk_variable_id),SAVE :: swr_id 52 PUBLIC trc_sms_fabm ! called by trcsms.F90 module 53 PUBLIC trc_sms_fabm_alloc ! called by trcini_fabm.F90 module 54 PUBLIC trc_sms_fabm_check_mass ! called by trcwri_fabm.F90 55 PUBLIC nemo_fabm_start 56 57 ! Work arrays 58 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: flux ! Cross-interface flux of pelagic variables (# m-2 s-1) 59 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: current_total ! Totals of conserved quantities 60 61 ! Arrays for environmental variables that are computed by the coupler 62 REAL(wp), PUBLIC, TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: prn,rho 63 REAL(wp), PUBLIC, TARGET, ALLOCATABLE, DIMENSION(:,:) :: taubot 64 REAL(wp), PUBLIC, TARGET :: daynumber_in_year 65 66 ! State repair counters 67 INTEGER, SAVE :: repair_interior_count = 0 68 INTEGER, SAVE :: repair_surface_count = 0 69 INTEGER, SAVE :: repair_bottom_count = 0 70 71 ! Coupler parameters 72 INTEGER, PUBLIC :: nn_adv ! Vertical advection scheme for sinking/floating/movement 73 ! (1: 1st order upwind, 3: 3rd order TVD) 74 75 ! Flag indicating whether model%start has been called (will be done on-demand) 76 LOGICAL, SAVE :: started = .false. 80 77 81 78 !!---------------------------------------------------------------------- … … 96 93 ! 97 94 INTEGER, INTENT(in) :: kt ! ocean time-step index 98 INTEGER :: jn 99 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrfabm 95 INTEGER :: jn, jk 96 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrfabm, pdat 97 REAL(wp), DIMENSION(jpi,jpj) :: vint 100 98 101 99 !!---------------------------------------------------------------------- … … 109 107 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 110 108 109 IF (.NOT. started) CALL nemo_fabm_start 110 111 111 CALL update_inputs( kt ) 112 112 113 CALL compute_fabm 114 115 CALL compute_vertical_movement( kt )113 CALL compute_fabm( kt ) 114 115 CALL compute_vertical_movement( kt, nn_adv ) 116 116 117 117 CALL st2d_fabm_nxt( kt ) … … 123 123 CALL trc_bc_read ( kt ) ! tracers: surface and lateral Boundary Conditions 124 124 CALL trc_rnf_fabm ( kt ) ! River forcings 125 126 ! Send 3D diagnostics to output (these apply to time "n") 127 DO jn = 1, size(model%interior_diagnostic_variables) 128 IF (model%interior_diagnostic_variables(jn)%save) THEN 129 ! Save 3D field 130 pdat => model%get_interior_diagnostic_data(jn) 131 CALL iom_put(model%interior_diagnostic_variables(jn)%name, pdat) 132 133 ! Save depth integral if selected for output in XIOS 134 IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT')) THEN 135 vint = 0._wp 136 DO jk = 1, jpkm1 137 vint = vint + pdat(:,:,jk) * fse3t(:,:,jk) * tmask(:,:,jk) 138 END DO 139 CALL iom_put(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT', vint) 140 END IF 141 END IF 142 END DO 143 144 ! Send 2D diagnostics to output (these apply to time "n") 145 DO jn = 1, size(model%horizontal_diagnostic_variables) 146 IF (model%horizontal_diagnostic_variables(jn)%save) & 147 CALL iom_put( model%horizontal_diagnostic_variables(jn)%name, model%get_horizontal_diagnostic_data(jn)) 148 END DO 125 149 126 150 IF( l_trdtrc ) THEN ! Save the trends in the mixed layer … … 139 163 INTEGER, INTENT(IN) :: kt 140 164 INTEGER :: ji,jj,jk,jkmax 141 REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d, zmld 165 REAL(wp), DIMENSION(jpi,jpj) :: zmld 166 REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d 142 167 143 168 IF (kt == nittrc000) THEN … … 148 173 PHYT_AVG(:,:) = 0.0 149 174 150 pgrow_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model,jp_fabm_pgrow)151 ploss_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model,jp_fabm_ploss)175 pgrow_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_pgrow) 176 ploss_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_ploss) 152 177 153 178 SELECT CASE( mld_choice_bgc ) … … 220 245 END SUBROUTINE asmdiags_fabm 221 246 222 SUBROUTINE compute_fabm() 247 SUBROUTINE compute_fabm( kt ) 248 INTEGER, INTENT(in) :: kt ! ocean time-step index 249 223 250 INTEGER :: ji,jj,jk,jn 224 TYPE(type_state) :: valid_state251 LOGICAL :: valid, repaired 225 252 REAL(wp) :: zalfg,zztmpx,zztmpy 226 253 227 254 ! Validate current model state (setting argument to .TRUE. enables repair=clipping) 228 valid_state = check_state(.TRUE.)229 IF (.NOT. valid_state%valid) THEN255 CALL check_state(.TRUE., valid, repaired) 256 IF (.NOT. valid) THEN 230 257 WRITE(numout,*) "Invalid value in FABM encountered in area ",narea,"!!!" 231 258 #if defined key_iomput … … 240 267 #endif 241 268 END IF 242 IF ( valid_state%repaired) THEN269 IF (repaired) THEN 243 270 WRITE(numout,*) "Total interior repairs up to now on process",narea,":",repair_interior_count 244 271 WRITE(numout,*) "Total surface repairs up to now on process",narea,":",repair_surface_count … … 246 273 ENDIF 247 274 275 daynumber_in_year = fjulday - fjulstartyear + 1 276 248 277 ! Compute the now hydrostatic pressure 249 278 ! copied from istate.F90 250 279 ! ------------------------------------ 251 280 252 zalfg = 0.5e-4 * grav ! FABM wants dbar, convert from Pa 253 254 rho = rau0 * ( 1. + rhd ) 255 256 prn(:,:,1) = 10.1325 + zalfg * fse3t(:,:,1) * rho(:,:,1) 257 258 daynumber_in_year=(fjulday-fjulstartyear+1)*1._wp 259 260 DO jk = 2, jpk ! Vertical integration from the surface 261 prn(:,:,jk) = prn(:,:,jk-1) + zalfg * ( & 262 fse3t(:,:,jk-1) * rho(:,:,jk-1) & 263 + fse3t(:,:,jk) * rho(:,:,jk) ) 264 END DO 265 266 ! Bottom stress 267 taubot(:,:) = 0._wp 268 DO jj = 2, jpjm1 269 DO ji = fs_2, fs_jpim1 ! vector opt. 270 zztmpx = ( bfrua(ji ,jj) * un(ji ,jj,mbku(ji ,jj)) & 271 & + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj)) ) 272 zztmpy = ( bfrva(ji, jj) * vn(ji,jj ,mbkv(ji,jj )) & 273 & + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1)) ) 274 taubot(ji,jj) = 0.5_wp * rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 275 ! 276 ENDDO 277 ENDDO 278 ! Compute light extinction 279 DO jk=1,jpk 280 DO jj=1,jpj 281 call fabm_get_light_extinction(model,1,jpi,jj,jk,ext) 282 END DO 283 END DO 284 285 ! Compute light field (stored among FABM's internal diagnostics) 286 DO jj=1,jpj 287 DO ji=1,jpi 288 call fabm_get_light(model,1,jpk,ji,jj) 289 END DO 290 END DO 291 292 ! TODO: retrieve 3D shortwave and store in etot3 281 IF (ALLOCATED(rho)) rho = rau0 * ( 1._wp + rhd ) 282 283 IF (ALLOCATED(prn)) THEN 284 zalfg = 0.5e-4_wp * grav ! FABM wants dbar, convert from Pa (and multiply with 0.5 to average 2 cell thicknesses below) 285 prn(:,:,1) = 10.1325_wp + zalfg * fse3t(:,:,1) * rho(:,:,1) 286 DO jk = 2, jpkm1 ! Vertical integration from the surface 287 prn(:,:,jk) = prn(:,:,jk-1) + zalfg * ( & 288 fse3t(:,:,jk-1) * rho(:,:,jk-1) & 289 + fse3t(:,:,jk) * rho(:,:,jk) ) 290 END DO 291 END IF 292 293 ! Compute the bottom stress 294 ! copied from diawri.F90 295 ! ------------------------------------ 296 297 IF (ALLOCATED(taubot)) THEN 298 taubot(:,:) = 0._wp 299 DO jj = 2, jpjm1 300 DO ji = fs_2, fs_jpim1 ! vector opt. 301 zztmpx = ( bfrua(ji ,jj) * un(ji ,jj,mbku(ji ,jj)) & 302 & + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj)) ) 303 zztmpy = ( bfrva(ji, jj) * vn(ji,jj ,mbkv(ji,jj )) & 304 & + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1)) ) 305 taubot(ji,jj) = 0.5_wp * rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 306 ! 307 END DO 308 END DO 309 END IF 310 311 CALL model%prepare_inputs(real(kt, wp),nyear,nmonth,nday,REAL(nsec_day,wp)) 312 313 ! Retrieve 3D shortwave and store in etot3 314 IF (ln_qsr_spec) THEN 315 etot3(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_swr) 316 ENDIF 293 317 294 318 ! Zero rate array of interface-attached state variables … … 296 320 297 321 ! Compute interfacial source terms and fluxes 298 DO jj= 1,jpj299 ! Process bottom ( fabm_do_bottomincrements rather than sets, so zero flux array first)322 DO jj=2,jpjm1 323 ! Process bottom (get_bottom_sources increments rather than sets, so zero flux array first) 300 324 flux = 0._wp 301 CALL fabm_do_bottom(model,1,jpi,jj,flux,fabm_st2Da(:,jj,jp_fabm_surface+1:))325 CALL model%get_bottom_sources(fs_2,fs_jpim1,jj,flux,fabm_st2Da(fs_2:fs_jpim1,jj,jp_fabm_surface+1:)) 302 326 DO jn=1,jp_fabm 303 DO ji=1,jpi 304 ! Divide bottom fluxes by height of bottom layer and add to source terms. 305 ! TODO: is there perhaps an existing variable for fse3t(ji,jj,mbkt(ji,jj))?? 306 tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) = tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) + flux(ji,jn)/fse3t(ji,jj,mbkt(ji,jj)) 307 END DO 308 END DO 309 310 ! Process surface (fabm_do_surface increments rather than sets, so zero flux array first) 327 ! Divide bottom fluxes by height of bottom layer and add to source terms. 328 DO ji=fs_2,fs_jpim1 329 tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) = tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) + flux(ji,jn)/fse3t(ji,jj,mbkt(ji,jj)) 330 END DO 331 END DO 332 333 ! Process surface (get_surface_sources increments rather than sets, so zero flux array first) 311 334 flux = 0._wp 312 CALL fabm_do_surface(model,1,jpi,jj,flux,fabm_st2Da(:,jj,1:jp_fabm_surface)) 335 CALL model%get_surface_sources(fs_2,fs_jpim1,jj,flux,fabm_st2Da(fs_2:fs_jpim1,jj,1:jp_fabm_surface)) 336 ! Divide surface fluxes by height of surface layer and add to source terms. 313 337 DO jn=1,jp_fabm 314 ! Divide surface fluxes by height of surface layer and add to source terms. 315 tra(:,jj,1,jp_fabm_m1+jn) = tra(:,jj,1,jp_fabm_m1+jn) + flux(:,jn)/fse3t(:,jj,1) 316 END DO 317 END DO 318 319 ! Compute interior source terms (NB fabm_do increments rather than sets) 320 DO jk=1,jpk 321 DO jj=1,jpj 322 CALL fabm_do(model,1,jpi,jj,jk,tra(:,jj,jk,jp_fabm0:jp_fabm1)) 323 END DO 324 END DO 338 DO ji=fs_2,fs_jpim1 339 tra(ji,jj,1,jp_fabm_m1+jn) = tra(ji,jj,1,jp_fabm_m1+jn) + flux(ji,jn)/fse3t(ji,jj,1) 340 END DO 341 END DO 342 END DO 343 344 ! Compute interior source terms (NB get_interior_sources increments rather than sets) 345 DO jk=1,jpkm1 346 DO jj=2,jpjm1 347 CALL model%get_interior_sources(fs_2,fs_jpim1,jj,jk,tra(fs_2:fs_jpim1,jj,jk,jp_fabm0:jp_fabm1)) 348 END DO 349 END DO 350 351 CALL model%finalize_outputs() 325 352 END SUBROUTINE compute_fabm 326 353 327 FUNCTION check_state(repair) RESULT(exit_state)328 LOGICAL, INTENT(IN) :: repair329 TYPE(type_state) :: exit_state330 331 INTEGER :: jj,jk332 LOGICAL :: valid_int,valid_sf,valid_bt333 334 exit_state%valid = .TRUE.335 exit_state%repaired =.FALSE.336 DO jk=1,jpk 337 DO jj= 1,jpj338 CALL fabm_check_state(model,1,jpi,jj,jk,repair,valid_int)339 IF (repair .AND..NOT.valid_int) THEN354 SUBROUTINE check_state(repair, valid, repaired) 355 LOGICAL, INTENT(IN) :: repair 356 LOGICAL, INTENT(OUT) :: valid, repaired 357 358 INTEGER :: jj, jk 359 LOGICAL :: valid_int, valid_sf, valid_bt 360 361 valid = .TRUE. ! Whether the model state is valid after this subroutine returns 362 repaired = .FALSE. ! Whether the model state has been repaired by this subroutine 363 DO jk=1,jpkm1 364 DO jj=2,jpjm1 365 CALL model%check_interior_state(fs_2, fs_jpim1, jj, jk, repair, valid_int) 366 IF (repair .AND. .NOT. valid_int) THEN 340 367 repair_interior_count = repair_interior_count + 1 341 exit_state%repaired = .TRUE.368 repaired = .TRUE. 342 369 END IF 343 IF (.NOT. (valid_int.OR.repair)) exit_state%valid = .FALSE.344 END DO 345 END DO 346 DO jj= 1,jpj347 CALL fabm_check_surface_state(model,1,jpi,jj,repair,valid_sf)348 IF (repair .AND..NOT.valid_sf) THEN370 IF (.NOT. (valid_int .OR. repair)) valid = .FALSE. 371 END DO 372 END DO 373 DO jj=2,jpjm1 374 CALL model%check_surface_state(fs_2, fs_jpim1, jj, repair, valid_sf) 375 IF (repair .AND. .NOT. valid_sf) THEN 349 376 repair_surface_count = repair_surface_count + 1 350 exit_state%repaired = .TRUE.377 repaired = .TRUE. 351 378 END IF 352 IF (.NOT. (valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE.353 CALL fabm_check_bottom_state(model,1,jpi,jj,repair,valid_bt)354 IF (repair .AND..NOT.valid_bt) THEN379 IF (.NOT. (valid_sf .AND. valid_bt) .AND. .NOT. repair) valid = .FALSE. 380 CALL model%check_bottom_state(fs_2, fs_jpim1, jj, repair, valid_bt) 381 IF (repair .AND. .NOT. valid_bt) THEN 355 382 repair_bottom_count = repair_bottom_count + 1 356 exit_state%repaired = .TRUE.383 repaired = .TRUE. 357 384 END IF 358 IF (.NOT. (valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE.359 END DO 360 END FUNCTION385 IF (.NOT. (valid_sf .AND. valid_bt) .AND. .NOT. repair) valid = .FALSE. 386 END DO 387 END SUBROUTINE 361 388 362 389 SUBROUTINE trc_sms_fabm_check_mass() 363 390 REAL(wp) :: total(SIZE(model%conserved_quantities)) 364 INTEGER :: j k,jj,jn391 INTEGER :: ji,jk,jj,jn 365 392 366 393 total = 0._wp 367 394 368 DO jk=1,jpk 369 DO jj=1,jpj 370 CALL fabm_get_conserved_quantities(model,1,jpi,jj,jk,current_total) 395 IF (.NOT. started) CALL nemo_fabm_start 396 397 DO jk=1,jpkm1 398 DO jj=2,jpjm1 399 CALL model%get_interior_conserved_quantities(fs_2,fs_jpim1,jj,jk,current_total) 371 400 DO jn=1,SIZE(model%conserved_quantities) 372 total(jn) = total(jn) + SUM(cvol(:,jj,jk)*current_total(:,jn)*tmask_i(:,jj)) 401 DO ji=fs_2,fs_jpim1 402 total(jn) = total(jn) + cvol(ji,jj,jk) * current_total(ji,jn) * tmask_i(ji,jj) 403 END DO 373 404 END DO 374 405 END DO 375 406 END DO 376 407 377 DO jj= 1,jpj378 CALL fabm_get_horizontal_conserved_quantities(model,1,jpi,jj,current_total)408 DO jj=2,jpjm1 409 CALL model%get_horizontal_conserved_quantities(fs_2,fs_jpim1,jj,current_total) 379 410 DO jn=1,SIZE(model%conserved_quantities) 380 total(jn) = total(jn) + SUM(e1e2t(:,jj)*current_total(:,jn)*tmask_i(:,jj)) 411 DO ji=fs_2,fs_jpim1 412 total(jn) = total(jn) + e1e2t(ji,jj) * current_total(ji,jn) * tmask_i(ji,jj) 413 END DO 381 414 END DO 382 415 END DO … … 434 467 435 468 INTEGER FUNCTION trc_sms_fabm_alloc() 436 INTEGER :: j j,jk,jn469 INTEGER :: jn 437 470 !!---------------------------------------------------------------------- 438 471 !! *** ROUTINE trc_sms_fabm_alloc *** … … 441 474 ! ALLOCATE here the arrays specific to FABM 442 475 ALLOCATE( lk_rad_fabm(jp_fabm)) 443 ALLOCATE(prn(jpi, jpj, jpk))444 ALLOCATE(rho(jpi, jpj, jpk))445 ALLOCATE(taubot(jpi, jpj))476 IF (model%variable_needs_values(fabm_standard_variables%pressure)) ALLOCATE(prn(jpi, jpj, jpk)) 477 IF (ALLOCATED(prn) .or. model%variable_needs_values(fabm_standard_variables%density)) ALLOCATE(rho(jpi, jpj, jpk)) 478 IF (model%variable_needs_values(fabm_standard_variables%bottom_stress)) ALLOCATE(taubot(jpi, jpj)) 446 479 ! ALLOCATE( tab(...) , STAT=trc_sms_fabm_alloc ) 447 480 … … 452 485 453 486 ! Work array to hold surface and bottom fluxes 454 ALLOCATE(flux(jpi,jp_fabm)) 455 456 ! Work array to hold extinction coefficients 457 ALLOCATE(ext(jpi)) 458 ext=0._wp 487 ALLOCATE(flux(fs_2:fs_jpim1,jp_fabm)) 459 488 460 489 ! Allocate work arrays for vertical movement 461 ALLOCATE(w_ct(jpi,jpk,jp_fabm)) 462 ALLOCATE(w_if(jpk,jp_fabm)) 463 ALLOCATE(zwgt_if(jpk,jp_fabm)) 464 ALLOCATE(flux_if(jpk,jp_fabm)) 465 ALLOCATE(current_total(jpi,SIZE(model%conserved_quantities))) 490 ALLOCATE(w_ct(fs_2:fs_jpim1,1:jpkm1,jp_fabm)) 491 ALLOCATE(current_total(fs_2:fs_jpim1,SIZE(model%conserved_quantities))) 466 492 #if defined key_trdtrc && defined key_iomput 467 493 IF( lk_trdtrc ) ALLOCATE(tr_vmv(jpi,jpj,jpk,jp_fabm)) … … 474 500 ! 475 501 476 ! Make FABM aware of diagnostics that are not needed [not included in output] 477 DO jn=1,size(model%diagnostic_variables) 478 !model%diagnostic_variables(jn)%save = iom_use(model%diagnostic_variables(jn)%name) 479 END DO 480 DO jn=1,size(model%horizontal_diagnostic_variables) 481 !model%horizontal_diagnostic_variables(jn)%save = iom_use(model%horizontal_diagnostic_variables(jn)%name) 482 END DO 483 484 ! Provide FABM with domain extents [after this, the save attribute of diagnostic variables can no longe change!] 485 call fabm_set_domain(model,jpi, jpj, jpk) 486 487 ! Provide FABM with the vertical indices of the surface and bottom, and the land-sea mask. 488 call model%set_bottom_index(mbkt) ! NB mbkt extents should match dimension lengths provided to fabm_set_domain 489 call model%set_surface_index(1) 490 call fabm_set_mask(model,tmask,tmask(:,:,1)) ! NB tmask extents should match dimension lengths provided to fabm_set_domain 491 492 ! Send pointers to state data to FABM 493 do jn=1,jp_fabm 494 call fabm_link_bulk_state_data(model,jn,trn(:,:,:,jp_fabm_m1+jn)) 495 end do 502 ! Provide FABM with domain extents 503 CALL model%set_domain(jpi, jpj, jpk, rdt) 504 CALL model%set_domain_start(fs_2, 2, 1) 505 CALL model%set_domain_stop(fs_jpim1, jpjm1, jpkm1) 506 507 ! Provide FABM with the vertical indices of the bottom, and the land-sea mask. 508 CALL model%set_bottom_index(mbkt) ! NB mbkt extents should match dimension lengths provided to set_domain 509 CALL model%set_mask(tmask,tmask(:,:,1)) ! NB tmask extents should match dimension lengths provided to set_domain 510 511 ! Initialize state and send pointers to state data to FABM 512 ! We mask land points in states with zeros, as per with NEMO "convention" 513 ! NB we cannot call model%initialize_*_state at this point, because model%start has not been called yet. 514 DO jn=1,jp_fabm 515 trn(:,:,:,jp_fabm_m1+jn) = model%interior_state_variables(jn)%initial_value * tmask 516 CALL model%link_interior_state_data(jn,trn(:,:,:,jp_fabm_m1+jn)) 517 END DO 496 518 DO jn=1,jp_fabm_surface 497 CALL fabm_link_surface_state_data(model,jn,fabm_st2Dn(:,:,jn)) 519 fabm_st2Dn(:,:,jn) = model%surface_state_variables(jn)%initial_value * tmask(:,:,1) 520 CALL model%link_surface_state_data(jn,fabm_st2Dn(:,:,jn)) 498 521 END DO 499 522 DO jn=1,jp_fabm_bottom 500 CALL fabm_link_bottom_state_data(model,jn,fabm_st2Dn(:,:,jp_fabm_surface+jn)) 523 fabm_st2Dn(:,:,jp_fabm_surface+jn) = model%bottom_state_variables(jn)%initial_value * tmask(:,:,1) 524 CALL model%link_bottom_state_data(jn,fabm_st2Dn(:,:,jp_fabm_surface+jn)) 501 525 END DO 502 526 503 527 ! Send pointers to environmental data to FABM 504 call fabm_link_bulk_data(model,standard_variables%temperature,tsn(:,:,:,jp_tem)) 505 call fabm_link_bulk_data(model,standard_variables%practical_salinity,tsn(:,:,:,jp_sal)) 506 call fabm_link_bulk_data(model,standard_variables%density,rho(:,:,:)) 507 call fabm_link_bulk_data(model,standard_variables%pressure,prn) 508 call fabm_link_horizontal_data(model,standard_variables%bottom_stress,taubot(:,:)) 509 ! correct target for cell thickness depends on NEMO configuration: 510 #ifdef key_vvl 511 call fabm_link_bulk_data(model,standard_variables%cell_thickness,e3t_n) 512 #else 513 call fabm_link_bulk_data(model,standard_variables%cell_thickness,e3t_0) 514 #endif 515 call fabm_link_horizontal_data(model,standard_variables%latitude,gphit) 516 call fabm_link_horizontal_data(model,standard_variables%longitude,glamt) 517 call fabm_link_scalar_data(model,standard_variables%number_of_days_since_start_of_the_year,daynumber_in_year) 518 call fabm_link_horizontal_data(model,standard_variables%wind_speed,wndm(:,:)) 519 call fabm_link_horizontal_data(model,standard_variables%surface_downwelling_shortwave_flux,qsr(:,:)) 520 call fabm_link_horizontal_data(model,standard_variables%bottom_depth_below_geoid,bathy(:,:)) 521 522 swr_id = model%get_bulk_variable_id(standard_variables%downwelling_shortwave_flux) 528 CALL model%link_interior_data(fabm_standard_variables%depth, fsdept(:,:,:)) 529 CALL model%link_interior_data(fabm_standard_variables%temperature, tsn(:,:,:,jp_tem)) 530 CALL model%link_interior_data(fabm_standard_variables%practical_salinity, tsn(:,:,:,jp_sal)) 531 IF (ALLOCATED(rho)) CALL model%link_interior_data(fabm_standard_variables%density, rho(:,:,:)) 532 IF (ALLOCATED(prn)) CALL model%link_interior_data(fabm_standard_variables%pressure, prn) 533 IF (ALLOCATED(taubot)) CALL model%link_horizontal_data(fabm_standard_variables%bottom_stress, taubot(:,:)) 534 CALL model%link_interior_data(fabm_standard_variables%cell_thickness, fse3t(:,:,:)) 535 CALL model%link_horizontal_data(fabm_standard_variables%latitude, gphit) 536 CALL model%link_horizontal_data(fabm_standard_variables%longitude, glamt) 537 CALL model%link_scalar(fabm_standard_variables%number_of_days_since_start_of_the_year, daynumber_in_year) 538 CALL model%link_horizontal_data(fabm_standard_variables%wind_speed, wndm(:,:)) 539 CALL model%link_horizontal_data(fabm_standard_variables%surface_downwelling_shortwave_flux, qsr(:,:)) 540 CALL model%link_horizontal_data(fabm_standard_variables%bottom_depth_below_geoid, bathy(:,:)) 541 CALL model%link_horizontal_data(fabm_standard_variables%ice_area_fraction, fr_i(:,:)) 523 542 524 543 ! Obtain user-specified input variables (read from NetCDF file) 525 call link_inputs 526 call update_inputs( nit000, .false. ) 527 528 ! Check whether FABM has all required data 529 call fabm_check_ready(model) 530 531 ! Initialize state 532 DO jj=1,jpj 533 CALL fabm_initialize_surface_state(model,1,jpi,jj) 534 CALL fabm_initialize_bottom_state(model,1,jpi,jj) 535 END DO 536 DO jk=1,jpk 537 DO jj=1,jpj 538 CALL fabm_initialize_state(model,1,jpi,jj,jk) 539 END DO 540 END DO 544 CALL link_inputs 545 CALL update_inputs(nit000, .FALSE.) 541 546 542 547 ! Set mask for negativity corrections to the relevant states 543 lk_rad_fabm = .FALSE.548 lk_rad_fabm(:) = .FALSE. 544 549 DO jn=1,jp_fabm 545 IF (model% state_variables(jn)%minimum.ge.0) THEN546 lk_rad_fabm(jn) =.TRUE.547 IF(lwp) WRITE(numout,*) 'FABM clipping for '//TRIM(model% state_variables(jn)%name)//' activated.'550 IF (model%interior_state_variables(jn)%minimum >= 0._wp) THEN 551 lk_rad_fabm(jn) = .TRUE. 552 IF(lwp) WRITE(numout,*) 'FABM clipping for '//TRIM(model%interior_state_variables(jn)%name)//' activated.' 548 553 END IF 549 554 END DO 550 555 551 ! Mask land points in states with zeros, not nice, but coherent552 ! with NEMO "convention":553 DO jn=jp_fabm0,jp_fabm1554 WHERE (tmask==0._wp)555 trn(:,:,:,jn)=0._wp556 END WHERE557 END DO558 DO jn=1,jp_fabm_surface+jp_fabm_bottom559 WHERE (tmask(:,:,1)==0._wp)560 fabm_st2Dn(:,:,jn)=0._wp561 END WHERE562 END DO563 556 564 557 ! Copy initial condition for interface-attached state variables to "previous" state field … … 566 559 fabm_st2Db = fabm_st2Dn 567 560 568 ! Initialise repair counters569 repair_interior_count = 0570 repair_surface_count = 0571 repair_bottom_count = 0572 573 561 END FUNCTION trc_sms_fabm_alloc 562 563 SUBROUTINE nemo_fabm_start() 564 INTEGER :: jn 565 566 ! Make FABM aware of diagnostics that are not needed [not included in output] 567 ! This works only after iom has completely initialised, because it depends on iom_use 568 DO jn=1,size(model%interior_diagnostic_variables) 569 model%interior_diagnostic_variables(jn)%save = iom_use(model%interior_diagnostic_variables(jn)%name) & 570 .or. iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT') & 571 .or. iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'25h') & 572 .or. iom_use('top_'//TRIM(model%interior_diagnostic_variables(jn)%name)) & 573 .or. iom_use('mid_'//TRIM(model%interior_diagnostic_variables(jn)%name)) & 574 .or. iom_use('bot_'//TRIM(model%interior_diagnostic_variables(jn)%name)) 575 END DO 576 model%interior_diagnostic_variables(jp_fabm_o3ta)%save = .TRUE. 577 model%interior_diagnostic_variables(jp_fabm_o3ph)%save = .TRUE. 578 model%interior_diagnostic_variables(jp_fabm_o3pc)%save = .TRUE. 579 model%interior_diagnostic_variables(jp_fabm_xeps)%save = .TRUE. 580 model%interior_diagnostic_variables(jp_fabm_pgrow)%save = .TRUE. 581 model%interior_diagnostic_variables(jp_fabm_ploss)%save = .TRUE. 582 IF( ln_qsr_spec ) THEN 583 model%interior_diagnostic_variables(jp_fabm_swr)%save = .TRUE. 584 ENDIF 585 DO jn=1,size(model%horizontal_diagnostic_variables) 586 model%horizontal_diagnostic_variables(jn)%save = iom_use(model%horizontal_diagnostic_variables(jn)%name) & 587 .or. iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h') 588 END DO 589 590 ! Check whether FABM has all required data 591 ! [after this, the save attribute of diagnostic variables can no longer change!] 592 CALL model%start() 593 594 started = .TRUE. 595 596 END SUBROUTINE 574 597 575 598 #else
Note: See TracChangeset
for help on using the changeset viewer.