Changeset 13241 for branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90
- Timestamp:
- 2020-07-03T14:42:49+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90
r10728 r13241 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 38 41 39 42 !USE fldread ! time interpolation 40 41 USE fabm42 43 43 44 IMPLICIT NONE … … 48 49 PRIVATE 49 50 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 51 PUBLIC trc_sms_fabm ! called by trcsms.F90 module 52 PUBLIC trc_sms_fabm_alloc ! called by trcini_fabm.F90 module 53 PUBLIC trc_sms_fabm_check_mass ! called by trcwri_fabm.F90 54 55 ! Work arrays 56 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: flux ! Cross-interface flux of pelagic variables (# m-2 s-1) 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: current_total ! Totals of conserved quantities 58 59 ! Arrays for environmental variables that are computed by the coupler 60 REAL(wp), PUBLIC, TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: prn,rho 61 REAL(wp), PUBLIC, TARGET, ALLOCATABLE, DIMENSION(:,:) :: taubot 62 REAL(wp), PUBLIC, TARGET :: daynumber_in_year 63 REAL(wp), PUBLIC, POINTER :: swrad(:,:,:) ! pointer to ERSEM spectral heating term 64 65 ! State repair counters 66 INTEGER, SAVE :: repair_interior_count = 0 67 INTEGER, SAVE :: repair_surface_count = 0 68 INTEGER, SAVE :: repair_bottom_count = 0 69 70 ! Coupler parameters 71 INTEGER, PUBLIC :: nn_adv ! Vertical advection scheme for sinking/floating/movement 72 ! (1: 1st order upwind, 3: 3rd order TVD) 73 74 ! Flag indicating whether model%start has been called (will be done on-demand) 75 LOGICAL, SAVE :: started = .false. 80 76 81 77 !!---------------------------------------------------------------------- … … 96 92 ! 97 93 INTEGER, INTENT(in) :: kt ! ocean time-step index 98 INTEGER :: jn 99 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrfabm 94 INTEGER :: jn, jk 95 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrfabm, pdat 96 REAL(wp), DIMENSION(jpi,jpj) :: vint 100 97 101 98 !!---------------------------------------------------------------------- … … 109 106 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 110 107 108 IF (.NOT. started) CALL nemo_fabm_start 109 111 110 CALL update_inputs( kt ) 112 111 113 CALL compute_fabm 114 115 CALL compute_vertical_movement( kt )112 CALL compute_fabm( kt ) 113 114 CALL compute_vertical_movement( kt, nn_adv ) 116 115 117 116 CALL st2d_fabm_nxt( kt ) … … 123 122 CALL trc_bc_read ( kt ) ! tracers: surface and lateral Boundary Conditions 124 123 CALL trc_rnf_fabm ( kt ) ! River forcings 124 125 ! Send 3D diagnostics to output (these apply to time "n") 126 DO jn = 1, size(model%interior_diagnostic_variables) 127 IF (model%interior_diagnostic_variables(jn)%save) THEN 128 ! Save 3D field 129 pdat => model%get_interior_diagnostic_data(jn) 130 CALL iom_put(model%interior_diagnostic_variables(jn)%name, pdat) 131 132 ! Save depth integral if selected for output in XIOS 133 IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT')) THEN 134 vint = 0._wp 135 DO jk = 1, jpkm1 136 vint = vint + pdat(:,:,jk) * fse3t(:,:,jk) * tmask(:,:,jk) 137 END DO 138 CALL iom_put(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT', vint) 139 END IF 140 END IF 141 END DO 142 143 ! Send 2D diagnostics to output (these apply to time "n") 144 DO jn = 1, size(model%horizontal_diagnostic_variables) 145 IF (model%horizontal_diagnostic_variables(jn)%save) & 146 CALL iom_put( model%horizontal_diagnostic_variables(jn)%name, model%get_horizontal_diagnostic_data(jn)) 147 END DO 125 148 126 149 IF( l_trdtrc ) THEN ! Save the trends in the mixed layer … … 139 162 INTEGER, INTENT(IN) :: kt 140 163 INTEGER :: ji,jj,jk,jkmax 141 REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d, zmld 164 REAL(wp), DIMENSION(jpi,jpj) :: zmld 165 REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d 142 166 143 167 IF (kt == nittrc000) THEN … … 148 172 PHYT_AVG(:,:) = 0.0 149 173 150 pgrow_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model,jp_fabm_pgrow)151 ploss_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model,jp_fabm_ploss)174 pgrow_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_pgrow) 175 ploss_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_ploss) 152 176 153 177 SELECT CASE( mld_choice_bgc ) … … 220 244 END SUBROUTINE asmdiags_fabm 221 245 222 SUBROUTINE compute_fabm() 246 SUBROUTINE compute_fabm( kt ) 247 INTEGER, INTENT(in) :: kt ! ocean time-step index 248 223 249 INTEGER :: ji,jj,jk,jn 224 TYPE(type_state) :: valid_state250 LOGICAL :: valid, repaired 225 251 REAL(wp) :: zalfg,zztmpx,zztmpy 226 252 227 253 ! Validate current model state (setting argument to .TRUE. enables repair=clipping) 228 valid_state = check_state(.TRUE.)229 IF (.NOT. valid_state%valid) THEN254 CALL check_state(.TRUE., valid, repaired) 255 IF (.NOT. valid) THEN 230 256 WRITE(numout,*) "Invalid value in FABM encountered in area ",narea,"!!!" 231 257 #if defined key_iomput … … 240 266 #endif 241 267 END IF 242 IF ( valid_state%repaired) THEN268 IF (repaired) THEN 243 269 WRITE(numout,*) "Total interior repairs up to now on process",narea,":",repair_interior_count 244 270 WRITE(numout,*) "Total surface repairs up to now on process",narea,":",repair_surface_count … … 246 272 ENDIF 247 273 274 daynumber_in_year = fjulday - fjulstartyear + 1 275 248 276 ! Compute the now hydrostatic pressure 249 277 ! copied from istate.F90 250 278 ! ------------------------------------ 251 279 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 280 IF (ALLOCATED(rho)) rho = rau0 * ( 1._wp + rhd ) 281 282 IF (ALLOCATED(prn)) THEN 283 zalfg = 0.5e-4_wp * grav ! FABM wants dbar, convert from Pa (and multiply with 0.5 to average 2 cell thicknesses below) 284 prn(:,:,1) = 10.1325_wp + zalfg * fse3t(:,:,1) * rho(:,:,1) 285 DO jk = 2, jpkm1 ! Vertical integration from the surface 286 prn(:,:,jk) = prn(:,:,jk-1) + zalfg * ( & 287 fse3t(:,:,jk-1) * rho(:,:,jk-1) & 288 + fse3t(:,:,jk) * rho(:,:,jk) ) 289 END DO 290 END IF 291 292 ! Compute the bottom stress 293 ! copied from diawri.F90 294 ! ------------------------------------ 295 296 IF (ALLOCATED(taubot)) THEN 297 taubot(:,:) = 0._wp 298 DO jj = 2, jpjm1 299 DO ji = fs_2, fs_jpim1 ! vector opt. 300 zztmpx = ( bfrua(ji ,jj) * un(ji ,jj,mbku(ji ,jj)) & 301 & + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj)) ) 302 zztmpy = ( bfrva(ji, jj) * vn(ji,jj ,mbkv(ji,jj )) & 303 & + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1)) ) 304 taubot(ji,jj) = 0.5_wp * rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 305 ! 306 END DO 307 END DO 308 END IF 309 310 CALL model%prepare_inputs(real(kt, wp),nyear,nmonth,nday,REAL(nsec_day,wp)) 311 312 ! Retrieve 3D shortwave and store in etot3 313 IF (ln_qsr_spec) etot3(:,:,:) = swrad(:,:,:) 293 314 294 315 ! Zero rate array of interface-attached state variables … … 296 317 297 318 ! Compute interfacial source terms and fluxes 298 DO jj= 1,jpj299 ! Process bottom ( fabm_do_bottomincrements rather than sets, so zero flux array first)319 DO jj=2,jpjm1 320 ! Process bottom (get_bottom_sources increments rather than sets, so zero flux array first) 300 321 flux = 0._wp 301 CALL fabm_do_bottom(model,1,jpi,jj,flux,fabm_st2Da(:,jj,jp_fabm_surface+1:))322 CALL model%get_bottom_sources(fs_2,fs_jpim1,jj,flux,fabm_st2Da(fs_2:fs_jpim1,jj,jp_fabm_surface+1:)) 302 323 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) 324 ! Divide bottom fluxes by height of bottom layer and add to source terms. 325 DO ji=fs_2,fs_jpim1 326 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)) 327 END DO 328 END DO 329 330 ! Process surface (get_surface_sources increments rather than sets, so zero flux array first) 311 331 flux = 0._wp 312 CALL fabm_do_surface(model,1,jpi,jj,flux,fabm_st2Da(:,jj,1:jp_fabm_surface)) 332 CALL model%get_surface_sources(fs_2,fs_jpim1,jj,flux,fabm_st2Da(fs_2:fs_jpim1,jj,1:jp_fabm_surface)) 333 ! Divide surface fluxes by height of surface layer and add to source terms. 313 334 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 335 DO ji=fs_2,fs_jpim1 336 tra(ji,jj,1,jp_fabm_m1+jn) = tra(ji,jj,1,jp_fabm_m1+jn) + flux(ji,jn)/fse3t(ji,jj,1) 337 END DO 338 END DO 339 END DO 340 341 ! Compute interior source terms (NB get_interior_sources increments rather than sets) 342 DO jk=1,jpkm1 343 DO jj=2,jpjm1 344 CALL model%get_interior_sources(fs_2,fs_jpim1,jj,jk,tra(fs_2:fs_jpim1,jj,jk,jp_fabm0:jp_fabm1)) 345 END DO 346 END DO 347 348 CALL model%finalize_outputs() 325 349 END SUBROUTINE compute_fabm 326 350 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) THEN351 SUBROUTINE check_state(repair, valid, repaired) 352 LOGICAL, INTENT(IN) :: repair 353 LOGICAL, INTENT(OUT) :: valid, repaired 354 355 INTEGER :: jj, jk 356 LOGICAL :: valid_int, valid_sf, valid_bt 357 358 valid = .TRUE. ! Whether the model state is valid after this subroutine returns 359 repaired = .FALSE. ! Whether the model state has been repaired by this subroutine 360 DO jk=1,jpkm1 361 DO jj=2,jpjm1 362 CALL model%check_interior_state(fs_2, fs_jpim1, jj, jk, repair, valid_int) 363 IF (repair .AND. .NOT. valid_int) THEN 340 364 repair_interior_count = repair_interior_count + 1 341 exit_state%repaired = .TRUE.365 repaired = .TRUE. 342 366 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) THEN367 IF (.NOT. (valid_int .OR. repair)) valid = .FALSE. 368 END DO 369 END DO 370 DO jj=2,jpjm1 371 CALL model%check_surface_state(fs_2, fs_jpim1, jj, repair, valid_sf) 372 IF (repair .AND. .NOT. valid_sf) THEN 349 373 repair_surface_count = repair_surface_count + 1 350 exit_state%repaired = .TRUE.374 repaired = .TRUE. 351 375 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) THEN376 IF (.NOT. (valid_sf .AND. valid_bt) .AND. .NOT. repair) valid = .FALSE. 377 CALL model%check_bottom_state(fs_2, fs_jpim1, jj, repair, valid_bt) 378 IF (repair .AND. .NOT. valid_bt) THEN 355 379 repair_bottom_count = repair_bottom_count + 1 356 exit_state%repaired = .TRUE.380 repaired = .TRUE. 357 381 END IF 358 IF (.NOT. (valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE.359 END DO 360 END FUNCTION382 IF (.NOT. (valid_sf .AND. valid_bt) .AND. .NOT. repair) valid = .FALSE. 383 END DO 384 END SUBROUTINE 361 385 362 386 SUBROUTINE trc_sms_fabm_check_mass() 363 387 REAL(wp) :: total(SIZE(model%conserved_quantities)) 364 INTEGER :: j k,jj,jn388 INTEGER :: ji,jk,jj,jn 365 389 366 390 total = 0._wp 367 391 368 DO jk=1,jpk 369 DO jj=1,jpj 370 CALL fabm_get_conserved_quantities(model,1,jpi,jj,jk,current_total) 392 IF (.NOT. started) CALL nemo_fabm_start 393 394 DO jk=1,jpkm1 395 DO jj=2,jpjm1 396 CALL model%get_interior_conserved_quantities(fs_2,fs_jpim1,jj,jk,current_total) 371 397 DO jn=1,SIZE(model%conserved_quantities) 372 total(jn) = total(jn) + SUM(cvol(:,jj,jk)*current_total(:,jn)*tmask_i(:,jj)) 373 END DO 374 END DO 375 END DO 376 377 DO jj=1,jpj 378 CALL fabm_get_horizontal_conserved_quantities(model,1,jpi,jj,current_total) 398 DO ji=fs_2,fs_jpim1 399 total(jn) = total(jn) + cvol(ji,jj,jk) * current_total(ji,jn) * tmask_i(ji,jj) 400 END DO 401 END DO 402 END DO 403 END DO 404 405 DO jj=2,jpjm1 406 CALL model%get_horizontal_conserved_quantities(fs_2,fs_jpim1,jj,current_total) 379 407 DO jn=1,SIZE(model%conserved_quantities) 380 total(jn) = total(jn) + SUM(e1e2t(:,jj)*current_total(:,jn)*tmask_i(:,jj)) 408 DO ji=fs_2,fs_jpim1 409 total(jn) = total(jn) + e1e2t(ji,jj) * current_total(ji,jn) * tmask_i(ji,jj) 410 END DO 381 411 END DO 382 412 END DO … … 434 464 435 465 INTEGER FUNCTION trc_sms_fabm_alloc() 436 INTEGER :: j j,jk,jn466 INTEGER :: jn 437 467 !!---------------------------------------------------------------------- 438 468 !! *** ROUTINE trc_sms_fabm_alloc *** … … 441 471 ! ALLOCATE here the arrays specific to FABM 442 472 ALLOCATE( lk_rad_fabm(jp_fabm)) 443 ALLOCATE(prn(jpi, jpj, jpk))444 ALLOCATE(rho(jpi, jpj, jpk))445 ALLOCATE(taubot(jpi, jpj))473 IF (model%variable_needs_values(fabm_standard_variables%pressure)) ALLOCATE(prn(jpi, jpj, jpk)) 474 IF (ALLOCATED(prn) .or. model%variable_needs_values(fabm_standard_variables%density)) ALLOCATE(rho(jpi, jpj, jpk)) 475 IF (model%variable_needs_values(fabm_standard_variables%bottom_stress)) ALLOCATE(taubot(jpi, jpj)) 446 476 ! ALLOCATE( tab(...) , STAT=trc_sms_fabm_alloc ) 447 477 … … 452 482 453 483 ! 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 484 ALLOCATE(flux(fs_2:fs_jpim1,jp_fabm)) 459 485 460 486 ! 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))) 487 ALLOCATE(w_ct(fs_2:fs_jpim1,1:jpkm1,jp_fabm)) 488 ALLOCATE(current_total(fs_2:fs_jpim1,SIZE(model%conserved_quantities))) 466 489 #if defined key_trdtrc && defined key_iomput 467 490 IF( lk_trdtrc ) ALLOCATE(tr_vmv(jpi,jpj,jpk,jp_fabm)) … … 474 497 ! 475 498 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 499 ! Provide FABM with domain extents 500 CALL model%set_domain(jpi, jpj, jpk, rdt) 501 CALL model%set_domain_start(fs_2, 2, 1) 502 CALL model%set_domain_stop(fs_jpim1, jpjm1, jpkm1) 503 504 ! Provide FABM with the vertical indices of the bottom, and the land-sea mask. 505 CALL model%set_bottom_index(mbkt) ! NB mbkt extents should match dimension lengths provided to set_domain 506 CALL model%set_mask(tmask,tmask(:,:,1)) ! NB tmask extents should match dimension lengths provided to set_domain 507 508 ! Initialize state and send pointers to state data to FABM 509 ! We mask land points in states with zeros, as per with NEMO "convention" 510 ! NB we cannot call model%initialize_*_state at this point, because model%start has not been called yet. 511 DO jn=1,jp_fabm 512 trn(:,:,:,jp_fabm_m1+jn) = model%interior_state_variables(jn)%initial_value * tmask 513 CALL model%link_interior_state_data(jn,trn(:,:,:,jp_fabm_m1+jn)) 514 END DO 496 515 DO jn=1,jp_fabm_surface 497 CALL fabm_link_surface_state_data(model,jn,fabm_st2Dn(:,:,jn)) 516 fabm_st2Dn(:,:,jn) = model%surface_state_variables(jn)%initial_value * tmask(:,:,1) 517 CALL model%link_surface_state_data(jn,fabm_st2Dn(:,:,jn)) 498 518 END DO 499 519 DO jn=1,jp_fabm_bottom 500 CALL fabm_link_bottom_state_data(model,jn,fabm_st2Dn(:,:,jp_fabm_surface+jn)) 520 fabm_st2Dn(:,:,jp_fabm_surface+jn) = model%bottom_state_variables(jn)%initial_value * tmask(:,:,1) 521 CALL model%link_bottom_state_data(jn,fabm_st2Dn(:,:,jp_fabm_surface+jn)) 501 522 END DO 502 523 503 524 ! 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) 525 CALL model%link_interior_data(fabm_standard_variables%depth, fsdept(:,:,:)) 526 CALL model%link_interior_data(fabm_standard_variables%temperature, tsn(:,:,:,jp_tem)) 527 CALL model%link_interior_data(fabm_standard_variables%practical_salinity, tsn(:,:,:,jp_sal)) 528 IF (ALLOCATED(rho)) CALL model%link_interior_data(fabm_standard_variables%density, rho(:,:,:)) 529 IF (ALLOCATED(prn)) CALL model%link_interior_data(fabm_standard_variables%pressure, prn) 530 IF (ALLOCATED(taubot)) CALL model%link_horizontal_data(fabm_standard_variables%bottom_stress, taubot(:,:)) 531 CALL model%link_interior_data(fabm_standard_variables%cell_thickness, fse3t(:,:,:)) 532 CALL model%link_horizontal_data(fabm_standard_variables%latitude, gphit) 533 CALL model%link_horizontal_data(fabm_standard_variables%longitude, glamt) 534 CALL model%link_scalar(fabm_standard_variables%number_of_days_since_start_of_the_year, daynumber_in_year) 535 CALL model%link_horizontal_data(fabm_standard_variables%wind_speed, wndm(:,:)) 536 CALL model%link_horizontal_data(fabm_standard_variables%surface_downwelling_shortwave_flux, qsr(:,:)) 537 CALL model%link_horizontal_data(fabm_standard_variables%bottom_depth_below_geoid, bathy(:,:)) 538 CALL model%link_horizontal_data(fabm_standard_variables%ice_area_fraction, fr_i(:,:)) 523 539 524 540 ! 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 541 CALL link_inputs 542 CALL update_inputs(nit000, .FALSE.) 541 543 542 544 ! Set mask for negativity corrections to the relevant states 543 lk_rad_fabm = .FALSE.545 lk_rad_fabm(:) = .FALSE. 544 546 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.'547 IF (model%interior_state_variables(jn)%minimum >= 0._wp) THEN 548 lk_rad_fabm(jn) = .TRUE. 549 IF(lwp) WRITE(numout,*) 'FABM clipping for '//TRIM(model%interior_state_variables(jn)%name)//' activated.' 548 550 END IF 549 551 END DO 550 552 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 553 564 554 ! Copy initial condition for interface-attached state variables to "previous" state field … … 566 556 fabm_st2Db = fabm_st2Dn 567 557 568 ! Initialise repair counters 569 repair_interior_count = 0 570 repair_surface_count = 0 571 repair_bottom_count = 0 558 ! Pointer to spectral heating term 559 swrad => model%get_data(model%get_interior_variable_id(standard_variables%net_rate_of_absorption_of_shortwave_energy_in_layer)) 560 572 561 573 562 END FUNCTION trc_sms_fabm_alloc 563 564 SUBROUTINE nemo_fabm_start() 565 INTEGER :: jn 566 567 ! Make FABM aware of diagnostics that are not needed [not included in output] 568 ! This works only after iom has completely initialised, because it depends on iom_use 569 DO jn=1,size(model%interior_diagnostic_variables) 570 model%interior_diagnostic_variables(jn)%save = iom_use(model%interior_diagnostic_variables(jn)%name) & 571 .or. iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT') 572 END DO 573 DO jn=1,size(model%horizontal_diagnostic_variables) 574 model%horizontal_diagnostic_variables(jn)%save = iom_use(model%horizontal_diagnostic_variables(jn)%name) 575 END DO 576 577 ! Check whether FABM has all required data 578 ! [after this, the save attribute of diagnostic variables can no longer change!] 579 CALL model%start() 580 581 started = .TRUE. 582 END SUBROUTINE 574 583 575 584 #else
Note: See TracChangeset
for help on using the changeset viewer.
