Changeset 13576
- Timestamp:
- 2020-10-09T12:35:11+02:00 (3 years ago)
- Location:
- branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO
- Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
r13318 r13576 1122 1122 #elif defined key_fabm 1123 1123 totalk_tavg(:,:,:) = totalk_tavg(:,:,:) + & 1124 & fabm_get_interior_diagnostic_data(model,jp_fabm_o3ta) / pnumtimes_tavg1124 & model%get_interior_diagnostic_data(jp_fabm_o3ta) / pnumtimes_tavg 1125 1125 totalk_tavg(:,:,:) = totalk_tavg(:,:,:) * tmask(:,:,:) 1126 1126 #endif … … 1167 1167 cchl_p_tavg(:,:,:) = cchl_p(:,:,:) 1168 1168 #elif defined key_fabm 1169 totalk_tavg(:,:,:) = fabm_get_interior_diagnostic_data(model,jp_fabm_o3ta)1169 totalk_tavg(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_o3ta) 1170 1170 totalk_tavg(:,:,:) = totalk_tavg(:,:,:) * tmask(:,:,:) 1171 1171 #endif -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
r10390 r13576 25 25 USE par_fabm 26 26 USE st2d_fabm, ONLY: fabm_st2dn 27 USE fabm, ONLY: fabm_get_interior_diagnostic_data, &28 & fabm_get_horizontal_diagnostic_data29 27 #endif 30 28 … … 211 209 END DO 212 210 DO jn = 1, jp_fabm_3d 213 fabm_3d_25h(:,:,:,jn) = fabm_get_interior_diagnostic_data(model, jn) 211 IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'25h')) THEN 212 fabm_3d_25h(:,:,:,jn) = model%get_interior_diagnostic_data(jn) 213 ENDIF 214 214 END DO 215 215 DO jn = 1, jp_fabm_surface … … 220 220 END DO 221 221 DO jn = 1, jp_fabm_2d 222 fabm_2d_25h(:,:,jn) = fabm_get_horizontal_diagnostic_data(model,jn) 222 IF (iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h')) THEN 223 fabm_2d_25h(:,:,jn) = model%get_horizontal_diagnostic_data(jn) 224 ENDIF 223 225 END DO 224 226 #endif … … 327 329 END DO 328 330 DO jn = 1, jp_fabm_3d 329 fabm_3d_25h(:,:,:,jn) = fabm_3d_25h(:,:,:,jn) + fabm_get_interior_diagnostic_data(model, jn) 331 IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'25h')) THEN 332 fabm_3d_25h(:,:,:,jn) = fabm_3d_25h(:,:,:,jn) + model%get_interior_diagnostic_data(jn) 333 ENDIF 330 334 END DO 331 335 DO jn = 1, jp_fabm_surface … … 336 340 END DO 337 341 DO jn = 1, jp_fabm_2d 338 fabm_2d_25h(:,:,jn) = fabm_2d_25h(:,:,jn) + fabm_get_horizontal_diagnostic_data(model,jn) 342 IF (iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h')) THEN 343 fabm_2d_25h(:,:,jn) = fabm_2d_25h(:,:,jn) + model%get_horizontal_diagnostic_data(jn) 344 ENDIF 339 345 END DO 340 346 #endif … … 401 407 DO jn = 1, jp_fabm 402 408 zw3d(:,:,:) = fabm_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 403 CALL iom_put( TRIM(model% state_variables(jn)%name)//"25h", zw3d )409 CALL iom_put( TRIM(model%interior_state_variables(jn)%name)//"25h", zw3d ) 404 410 END DO 405 411 DO jn = 1, jp_fabm_3d 406 zw3d(:,:,:) = fabm_3d_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 407 CALL iom_put( TRIM(model%diagnostic_variables(jn)%name)//"25h", zw3d ) 412 IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'25h')) THEN 413 zw3d(:,:,:) = fabm_3d_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 414 CALL iom_put( TRIM(model%interior_diagnostic_variables(jn)%name)//"25h", zw3d ) 415 ENDIF 408 416 END DO 409 417 DO jn = 1, jp_fabm_surface … … 416 424 END DO 417 425 DO jn = 1, jp_fabm_2d 418 zw2d(:,:) = fabm_2d_25h(:,:,jn)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 419 CALL iom_put( TRIM(model%horizontal_diagnostic_variables(jn)%name)//"25h", zw2d ) 426 IF (iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h')) THEN 427 zw2d(:,:) = fabm_2d_25h(:,:,jn)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 428 CALL iom_put( TRIM(model%horizontal_diagnostic_variables(jn)%name)//"25h", zw2d ) 429 ENDIF 420 430 END DO 421 431 #endif … … 468 478 END DO 469 479 DO jn = 1, jp_fabm_3d 470 fabm_3d_25h(:,:,:,jn) = fabm_get_interior_diagnostic_data(model, jn) 480 IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'25h')) THEN 481 fabm_3d_25h(:,:,:,jn) = model%get_interior_diagnostic_data(jn) 482 ENDIF 471 483 END DO 472 484 DO jn = 1, jp_fabm_surface … … 477 489 END DO 478 490 DO jn = 1, jp_fabm_2d 479 fabm_2d_25h(:,:,jn) = fabm_get_horizontal_diagnostic_data(model,jn) 491 IF (iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h')) THEN 492 fabm_2d_25h(:,:,jn) = model%get_horizontal_diagnostic_data(jn) 493 ENDIF 480 494 END DO 481 495 #endif -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90
r10390 r13576 14 14 USE trc, ONLY: trn 15 15 USE par_fabm 16 USE fabm, ONLY: fabm_get_interior_diagnostic_data17 16 #endif 18 17 … … 172 171 DO jn = 1, jp_fabm 173 172 CALL dia_calctmb( trn(:,:,:,jp_fabm_m1+jn), zwtmb ) 174 CALL iom_put( "top_"//TRIM(model% state_variables(jn)%name) , zwtmb(:,:,1) )175 CALL iom_put( "mid_"//TRIM(model% state_variables(jn)%name) , zwtmb(:,:,2) )176 CALL iom_put( "bot_"//TRIM(model% state_variables(jn)%name) , zwtmb(:,:,3) )173 CALL iom_put( "top_"//TRIM(model%interior_state_variables(jn)%name) , zwtmb(:,:,1) ) 174 CALL iom_put( "mid_"//TRIM(model%interior_state_variables(jn)%name) , zwtmb(:,:,2) ) 175 CALL iom_put( "bot_"//TRIM(model%interior_state_variables(jn)%name) , zwtmb(:,:,3) ) 177 176 END DO 178 177 DO jn = 1, jp_fabm_3d 179 CALL dia_calctmb( fabm_get_interior_diagnostic_data(model, jn), zwtmb ) 180 CALL iom_put( "top_"//TRIM(model%diagnostic_variables(jn)%name) , zwtmb(:,:,1) ) 181 CALL iom_put( "mid_"//TRIM(model%diagnostic_variables(jn)%name) , zwtmb(:,:,2) ) 182 CALL iom_put( "bot_"//TRIM(model%diagnostic_variables(jn)%name) , zwtmb(:,:,3) ) 178 IF ( iom_use('top_'//TRIM(model%interior_diagnostic_variables(jn)%name)) .OR. & 179 & iom_use('mid_'//TRIM(model%interior_diagnostic_variables(jn)%name)) .OR. & 180 & iom_use('bot_'//TRIM(model%interior_diagnostic_variables(jn)%name)) ) THEN 181 CALL dia_calctmb( model%get_interior_diagnostic_data(jn), zwtmb ) 182 CALL iom_put( "top_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,1) ) 183 CALL iom_put( "mid_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,2) ) 184 CALL iom_put( "bot_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,3) ) 185 ENDIF 183 186 END DO 184 187 #endif -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r8059 r13576 32 32 USE wrk_nemo ! Memory Allocation 33 33 USE timing ! Timing 34 #if defined key_fabm 35 USE trc, ONLY: trn ! FABM variables 36 USE par_fabm ! FABM parameters 37 #endif 34 38 35 39 IMPLICIT NONE … … 45 49 LOGICAL , PUBLIC :: ln_qsr_bio !: bio-model light absorption flag 46 50 LOGICAL , PUBLIC :: ln_qsr_ice !: light penetration for ice-model LIM3 (clem) 51 LOGICAL , PUBLIC :: ln_qsr_spec !: spectral model heating from ERSEM 47 52 INTEGER , PUBLIC :: nn_chldta !: use Chlorophyll data (=1) or not (=0) 48 53 INTEGER , PUBLIC :: nn_kd490dta !: use kd490dta data (=1) or not (=0) … … 106 111 REAL(wp) :: zz0, zz1, z1_e3t ! - - 107 112 REAL(wp), POINTER, DIMENSION(:,: ) :: zekb, zekg, zekr 113 REAL(wp), POINTER, DIMENSION(:,:,:) :: zekb_3d, zekg_3d, zekr_3d 108 114 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt 109 115 !!---------------------------------------------------------------------- … … 113 119 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 114 120 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 121 CALL wrk_alloc( jpi, jpj, jpk, zekb_3d, zekg_3d, zekr_3d ) 115 122 ! 116 123 IF( kt == nit000 ) THEN … … 150 157 151 158 ! ! ============================================== ! 152 IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! 159 IF( ln_qsr_spec ) THEN ! ERSEM spectral heating ! 160 ! ! ============================================== ! 161 DO jk = 1, jpkm1 162 qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) ) 163 END DO 164 ! Add to the general trend 165 DO jk = 1, jpkm1 166 DO jj = 2, jpjm1 167 DO ji = fs_2, fs_jpim1 ! vector opt. 168 z1_e3t = zfact / fse3t(ji,jj,jk) 169 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t 170 END DO 171 END DO 172 END DO 173 zea(:,:,:) = etot3(:,:,:) * tmask(:,:,:) 174 DO jk = jpkm1, 1, -1 175 zea(:,:,jk) = zea(:,:,jk) + zea(:,:,jk+1) 176 END DO 177 CALL iom_put( 'qsr3d', zea ) ! Shortwave Radiation 3D distribution 178 IF ( ln_qsr_ice ) THEN 179 DO jj = 1, jpj 180 DO ji = 1, jpi 181 IF ( qsr(ji,jj) /= 0._wp ) THEN 182 fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) ) 183 ELSE 184 fraqsr_1lev(ji,jj) = 1. 185 ENDIF 186 END DO 187 END DO 188 ENDIF 189 ! 190 191 ! ! ============================================== ! 192 ELSEIF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! 153 193 ! ! ============================================== ! 154 194 DO jk = 1, jpkm1 … … 185 225 ! ! ------------------------- ! 186 226 ! Set chlorophyl concentration 187 IF( nn_chldta == 1 .OR. lk_vvl ) THEN !* Variable Chlorophyll or ocean volume 188 ! 189 IF( nn_chldta == 1 ) THEN !* Variable Chlorophyll 227 IF( nn_chldta == 2 .OR. nn_chldta == 1 .OR. lk_vvl ) THEN !* Variable Chlorophyll or ocean volume 228 ! 229 IF( nn_chldta == 2 ) THEN 230 DO jk = 1, nksr+1 231 DO jj = 1, jpj 232 DO ji = 1, jpi 233 #if defined key_fabm 234 zchl = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) + & 235 & trn(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) 236 #endif 237 zchl = MIN( 10. , MAX( 0.03, zchl ) ) 238 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 239 ! 240 zekb_3d(ji,jj,jk) = rkrgb(1,irgb) 241 zekg_3d(ji,jj,jk) = rkrgb(2,irgb) 242 zekr_3d(ji,jj,jk) = rkrgb(3,irgb) 243 END DO 244 END DO 245 END DO 246 ! 247 ELSEIF( nn_chldta == 1 ) THEN !* Variable Chlorophyll 190 248 ! 191 249 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step … … 219 277 ! 220 278 DO jk = 2, nksr+1 279 IF( nn_chldta == 2 ) THEN 280 zekb(:,:) = zekb_3d(:,:,jk) 281 zekg(:,:) = zekg_3d(:,:,jk) 282 zekr(:,:) = zekr_3d(:,:,jk) 283 ENDIF 221 284 !CDIR NOVERRCHK 222 285 DO jj = 1, jpj … … 237 300 ! clem: store attenuation coefficient of the first ocean level 238 301 IF ( ln_qsr_ice ) THEN 302 IF( nn_chldta == 2 ) THEN 303 zekb(:,:) = zekb_3d(:,:,1) 304 zekg(:,:) = zekg_3d(:,:,1) 305 zekr(:,:) = zekr_3d(:,:,1) 306 ENDIF 239 307 DO jj = 1, jpj 240 308 DO ji = 1, jpi … … 386 454 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 387 455 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 456 CALL wrk_dealloc( jpi, jpj, jpk, zekb_3d, zekg_3d, zekr_3d ) 388 457 ! 389 458 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr') … … 423 492 !! 424 493 NAMELIST/namtra_qsr/ sn_chl, sn_kd490, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice, & 425 & nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta494 & ln_qsr_spec, nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta 426 495 !!---------------------------------------------------------------------- 427 496 … … 451 520 WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd 452 521 WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio 522 WRITE(numout,*) ' ERSEM spectral heating model ln_qsr_spec= ', ln_qsr_spec 453 523 WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice 454 WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta= ', nn_chldta524 WRITE(numout,*) ' RGB: model (2), file (1) or cst (0) chl nn_chldta = ', nn_chldta 455 525 WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs 456 526 WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0 … … 470 540 IF( ln_qsr_2bd ) ioptio = ioptio + 1 471 541 IF( ln_qsr_bio ) ioptio = ioptio + 1 542 IF( ln_qsr_spec ) ioptio = ioptio + 1 472 543 IF( nn_kd490dta == 1 ) ioptio = ioptio + 1 473 544 ! … … 478 549 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1 479 550 IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = 2 480 IF( ln_qsr_2bd ) nqsr = 3 481 IF( ln_qsr_bio ) nqsr = 4 482 IF( nn_kd490dta == 1 ) nqsr = 5 551 IF( ln_qsr_rgb .AND. nn_chldta == 2 ) nqsr = 3 552 IF( ln_qsr_2bd ) nqsr = 4 553 IF( ln_qsr_bio ) nqsr = 5 554 IF( nn_kd490dta == 1 ) nqsr = 6 555 IF( ln_qsr_spec ) nqsr = 7 483 556 ! 484 557 IF(lwp) THEN ! Print the choice 485 558 WRITE(numout,*) 486 559 IF( nqsr == 1 ) WRITE(numout,*) ' R-G-B light penetration - Constant Chlorophyll' 487 IF( nqsr == 2 ) WRITE(numout,*) ' R-G-B light penetration - Chl data ' 488 IF( nqsr == 3 ) WRITE(numout,*) ' 2 bands light penetration' 489 IF( nqsr == 4 ) WRITE(numout,*) ' bio-model light penetration' 490 IF( nqsr == 5 ) WRITE(numout,*) ' KD490 light penetration' 491 ENDIF 560 IF( nqsr == 2 ) WRITE(numout,*) ' R-G-B light penetration - Chl data from file' 561 IF( nqsr == 3 ) WRITE(numout,*) ' R-G-B light penetration - Chl data from model' 562 IF( nqsr == 4 ) WRITE(numout,*) ' 2 bands light penetration' 563 IF( nqsr == 5 ) WRITE(numout,*) ' bio-model light penetration' 564 IF( nqsr == 6 ) WRITE(numout,*) ' KD490 light penetration' 565 IF( nqsr == 7 ) WRITE(numout,*) ' ERSEM spectral light penetration' 566 ENDIF 567 #if ! defined key_fabm 568 ! 569 IF( nqsr == 2 ) THEN 570 CALL ctl_stop( 'nn_chldta=2 so trying to use ERSEM chlorophyll for light penetration', & 571 & 'but not running with ERSEM' ) 572 ELSEIF( nqsr == 7 ) THEN 573 CALL ctl_stop( 'ln_qsr_spec=.true. so trying to use ERSEM spectral light penetration', & 574 & 'but not running with ERSEM' ) 575 ENDIF 576 #endif 492 577 ! 493 578 ENDIF … … 536 621 & 'Solar penetration function of read chlorophyll', 'namtra_qsr' ) 537 622 ! 623 ELSEIF( nn_chldta == 2 ) THEN !* Chl data will be got from model at each time step 624 IF(lwp) WRITE(numout,*) 625 IF(lwp) WRITE(numout,*) ' Chlorophyll will be taken from model at each time step' 538 626 ELSE !* constant Chl : compute once for all the distribution of light (etot3) 539 627 IF(lwp) WRITE(numout,*) -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r11277 r13576 158 158 END DO 159 159 #else 160 IF( lk_asminc ) THEN 161 IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 ) ! Output background fields 162 IF( ln_asmdin ) THEN ! Direct initialization 163 IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 ) ! Tracers 164 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 ) ! Dynamics 165 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 ) ! SSH 166 IF( lk_bgcinc ) CALL bgc_asm_inc( nit000 - 1 ) ! BGC 167 ENDIF 168 ENDIF 160 ! Initial call to assimilation routines moved to stp 169 161 170 162 #if defined key_agrif … … 493 485 IF( lk_diaobs ) THEN ! Observation & model comparison 494 486 CALL dia_obs_init ! Initialize observational data 495 CALL dia_obs( nit000 - 1 ) ! Observation operator for restart 496 ENDIF 497 498 ! ! Assimilation increments 499 IF( lk_asminc ) THEN 500 #if defined key_shelf 501 CALL zdf_mxl_tref() ! Initialization of hmld_tref 502 #endif 503 CALL asm_inc_init ! Initialize assimilation increments 504 ENDIF 487 ! Initial call to dia_obs moved to stp 488 ENDIF 489 490 ! Initialisation of assimilation moved to stp 505 491 506 492 IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler 507 CALL dia_tmb_init ! TMB outputs508 CALL dia_25h_init ! 25h mean outputs509 CALL dia_diaopfoam_init ! FOAM operational output493 494 ! Initialisation of tmb/25h/diaopfoam outputs moved to stp 495 510 496 ! 511 497 END SUBROUTINE nemo_init -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/step.F90
r11277 r13576 96 96 CALL iom_init( cxios_context ) ! iom_put initialization 97 97 IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! initialize context for coarse grid 98 99 #if defined key_fabm 100 ! FABM can only finish initialising once IOM has 101 IF ( lk_fabm ) CALL nemo_fabm_start 102 #endif 103 104 ! First call to dia_* and asm_inc_init must wait for FABM to be initialised 105 ! (if running with FABM, but no harm moving the calls to here from nemo_init either way) 106 CALL dia_obs( nit000 - 1 ) ! Observation operator for restart 107 108 IF( lk_asminc ) THEN 109 #if defined key_shelf 110 CALL zdf_mxl_tref() ! Initialization of hmld_tref 111 #endif 112 CALL asm_inc_init ! Initialize assimilation increments 113 IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 ) ! Output background fields 114 IF( ln_asmdin ) THEN ! Direct initialization 115 IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 ) ! Tracers 116 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 ) ! Dynamics 117 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 ) ! SSH 118 IF( lk_bgcinc ) CALL bgc_asm_inc( nit000 - 1 ) ! BGC 119 ENDIF 120 ENDIF 121 122 CALL dia_tmb_init ! TMB outputs 123 CALL dia_25h_init ! 25h mean outputs 124 CALL dia_diaopfoam_init ! FOAM operational output 98 125 ENDIF 99 126 -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r11277 r13576 126 126 USE trcstp ! passive tracer time-stepping (trc_stp routine) 127 127 #endif 128 #if defined key_fabm 129 USE par_fabm, ONLY: & ! FABM parameters 130 & lk_fabm 131 USE trcsms_fabm, ONLY: & ! FABM routines 132 & nemo_fabm_start 133 #endif 134 USE diatmb ! Top,middle,bottom output 135 USE dia25h ! 25h mean output 136 USE diaopfoam ! FOAM operational output 128 137 !!---------------------------------------------------------------------- 129 138 !! NEMO/OPA 3.3 , NEMO Consortium (2010) -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90
r10158 r13576 19 19 USE fldread 20 20 USE par_fabm 21 USE fabm 21 USE fabm, only: type_fabm_horizontal_variable_id 22 22 23 23 IMPLICIT NONE 24 25 # include "vectopt_loop_substitute.h90" 24 26 25 27 PRIVATE … … 40 42 41 43 TYPE, PUBLIC, EXTENDS(type_input_variable) :: type_input_data 42 TYPE(type_ horizontal_variable_id):: horizontal_id43 TYPE(type_input_data), POINTER :: next => null()44 TYPE(type_fabm_horizontal_variable_id) :: horizontal_id 45 TYPE(type_input_data), POINTER :: next => null() 44 46 END TYPE 45 47 TYPE (type_input_data), POINTER, PUBLIC :: first_input_data => NULL() … … 87 89 ALLOCATE(input_data, STAT=ierr) 88 90 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'inputs_fabm:initialize_inputs: unable to allocate input_data object for variable '//TRIM(name) ) 89 input_data%horizontal_id = fabm_get_horizontal_variable_id(model,name)90 IF (.NOT. fabm_is_variable_used(input_data%horizontal_id)) THEN91 input_data%horizontal_id = model%get_horizontal_variable_id(name) 92 IF (.NOT.model%is_variable_used(input_data%horizontal_id)) THEN 91 93 ! This variable was not found among FABM's horizontal variables (at least, those that are read by one or more FABM modules) 92 94 CALL ctl_stop('STOP', 'inputs_fabm:initialize_inputs: variable "'//TRIM(name)//'" was not found among horizontal FABM variables.') … … 130 132 ! within tracer field 131 133 DO jn=1,jp_fabm 132 IF (TRIM(name) == TRIM(model% state_variables(jn)%name)) THEN134 IF (TRIM(name) == TRIM(model%interior_state_variables(jn)%name)) THEN 133 135 river_data%jp_pos = jp_fabm_m1+jn 134 136 END IF … … 173 175 ! Provide FABM with pointer to field that will receive prescribed data. 174 176 ! NB source=data_source_user guarantees that the prescribed data takes priority over any data FABM may already have for that variable. 175 CALL fabm_link_horizontal_data(model,input_data%horizontal_id,input_data%sf(1)%fnow(:,:,1),source=data_source_user)177 CALL model%link_horizontal_data(input_data%horizontal_id,input_data%sf(1)%fnow(:,:,1),source=data_source_user) 176 178 input_data => input_data%next 177 179 END DO … … 226 228 #endif 227 229 IF( kt == nit000 .OR. ( kt /= nit000 ) ) THEN 228 DO jj = 1, jpj229 DO ji = 1, jpi230 DO jj = 2, jpjm1 231 DO ji = fs_2, fs_jpim1 230 232 ! convert units and divide by surface area 231 233 ! loading / cell volume * vertical fraction of riverload -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90
r11545 r13576 1 1 MODULE par_fabm 2 2 3 #if defined key_fabm 4 # include "fabm_version.h" 5 # if _FABM_API_VERSION_ < 1 6 # error You need FABM 1.0 or later 7 # endif 3 8 USE fabm 9 #endif 4 10 5 11 IMPLICIT NONE 6 7 TYPE (type_model) :: model !FABM model instance8 12 9 13 INTEGER, PUBLIC :: jp_fabm0, jp_fabm1, jp_fabm, & … … 35 39 jp_fabm_r8c, jp_fabm_r8p, & 36 40 jp_fabm_r8s, jp_fabm_xeps, & 41 jp_fabm_swr, jp_fabm_kd490, & 37 42 jp_fabm_pgrow, jp_fabm_ploss 38 43 44 LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) :: lk_rad_fabm !: FABM negativity correction flag array 45 39 46 #if defined key_fabm 47 CLASS (type_fabm_model), POINTER :: model !FABM model instance 48 40 49 !!--------------------------------------------------------------------- 41 50 !! 'key_fabm' FABM tracers 42 51 !!--------------------------------------------------------------------- 43 52 LOGICAL, PUBLIC, PARAMETER :: lk_fabm = .TRUE. !: FABM flag 44 LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) :: lk_rad_fabm !: FABM negativity correction flag array45 53 #else 46 54 !!--------------------------------------------------------------------- -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90
r12352 r13576 4 4 !! TOP : initialisation of the FABM tracers 5 5 !!====================================================================== 6 !! History : 2.0 ! 2007-12 (C. Ethe, G. Madec) Original code 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 … … 17 18 USE par_fabm 18 19 USE trcsms_fabm 19 USE fabm _config,ONLY: fabm_create_model_from_yaml_file20 USE fabm ,ONLY: type_external_variable, fabm_initialize_library21 USE inputs_fabm,ONLY: initialize_inputs, link_inputs, &20 USE fabm, only: fabm_create_model, type_fabm_variable 21 USE fabm_driver 22 USE inputs_fabm,ONLY: initialize_inputs, link_inputs, & 22 23 type_input_variable,type_input_data,type_river_data, & 23 24 first_input_data,first_river_data … … 33 34 34 35 #if defined key_git_version 35 !#include "gitversion.h90"36 # include "gitversion.h90" 36 37 CHARACTER(len=*),parameter :: git_commit_id = _NEMO_COMMIT_ID_ 37 38 CHARACTER(len=*),parameter :: git_branch_name = _NEMO_BRANCH_ … … 39 40 40 41 PUBLIC trc_ini_fabm ! called by trcini.F90 module 41 PUBLIC nemo_fabm_init 42 PUBLIC nemo_fabm_configure 43 44 TYPE,extends(type_base_driver) :: type_nemo_fabm_driver 45 contains 46 procedure :: fatal_error => nemo_fabm_driver_fatal_error 47 procedure :: log_message => nemo_fabm_driver_log_message 48 end type 42 49 43 50 !!---------------------------------------------------------------------- … … 48 55 CONTAINS 49 56 50 SUBROUTINE nemo_fabm_ init()57 SUBROUTINE nemo_fabm_configure() 51 58 INTEGER :: jn 52 59 INTEGER, PARAMETER :: xml_unit = 1979 … … 55 62 CLASS (type_input_variable),POINTER :: input_pointer 56 63 64 ALLOCATE(type_nemo_fabm_driver::driver) 65 57 66 ! Allow FABM to parse fabm.yaml. This ensures numbers of variables are known. 58 call fabm_create_model_from_yaml_file(model)59 60 jp_fabm = size(model% state_variables)67 model => fabm_create_model() 68 69 jp_fabm = size(model%interior_state_variables) 61 70 jp_fabm_bottom = size(model%bottom_state_variables) 62 71 jp_fabm_surface = size(model%surface_state_variables) … … 66 75 jptra = jptra + jp_fabm 67 76 jp_fabm_2d = size(model%horizontal_diagnostic_variables) 68 jp_fabm_3d = size(model% diagnostic_variables)77 jp_fabm_3d = size(model%interior_diagnostic_variables) 69 78 jpdia2d = jpdia2d + jp_fabm_2d 70 79 jpdia3d = jpdia3d + jp_fabm_3d 71 80 jpdiabio = jpdiabio + jp_fabm 72 81 73 !Initialize input data structures. 82 ! Read inputs (river and additional 2D forcing) from fabm_input.nml 83 ! This must be done before writing field_def_fabm.xml, as that file 84 ! also describes the additional input variables. 74 85 call initialize_inputs 75 86 … … 123 134 jp_fabm_o3pc = fabm_diag_index( 'O3_pCO2' ) 124 135 jp_fabm_xeps = fabm_diag_index( 'light_xEPS' ) 136 jp_fabm_swr = fabm_diag_index( 'light_swr_abs' ) 137 jp_fabm_kd490 = fabm_diag_index( 'light_Kd_band3' ) 125 138 jp_fabm_pgrow = fabm_diag_index( 'p_grow_sum_result' ) 126 139 jp_fabm_ploss = fabm_diag_index( 'p_loss_sum_result' ) … … 128 141 IF (lwp) THEN 129 142 ! write field_def_fabm.xml on lead process 130 OPEN(UNIT=xml_unit, FILE='field_def_fabm.xml',ACTION='WRITE',STATUS='REPLACE')143 OPEN(UNIT=xml_unit, FILE='field_def_fabm.xml', ACTION='WRITE', STATUS='REPLACE') 131 144 132 145 WRITE (xml_unit,1000) '<field_definition level="1" prec="4" operation="average" enabled=".TRUE." default_value="1.e20" >' … … 134 147 WRITE (xml_unit,1000) ' <field_group id="ptrc_T" grid_ref="grid_T_3D">' 135 148 DO jn=1,jp_fabm 136 CALL write_variable_xml(xml_unit,model% state_variables(jn))149 CALL write_variable_xml(xml_unit,model%interior_state_variables(jn)) 137 150 #if defined key_trdtrc 138 CALL write_trends_xml(xml_unit,model% state_variables(jn))139 #endif 140 CALL write_25hourm_xml(xml_unit,model% state_variables(jn))141 CALL write_tmb_xml(xml_unit,model% state_variables(jn))151 CALL write_trends_xml(xml_unit,model%interior_state_variables(jn)) 152 #endif 153 CALL write_25hourm_xml(xml_unit,model%interior_state_variables(jn)) 154 CALL write_tmb_xml(xml_unit,model%interior_state_variables(jn)) 142 155 END DO 143 156 WRITE (xml_unit,1000) ' </field_group>' … … 155 168 156 169 WRITE (xml_unit,1000) ' <field_group id="diad_T" grid_ref="grid_T_2D">' 157 DO jn=1,size(model% diagnostic_variables)158 CALL write_variable_xml(xml_unit,model% diagnostic_variables(jn),3)159 CALL write_25hourm_xml(xml_unit,model% diagnostic_variables(jn),3)160 CALL write_tmb_xml(xml_unit,model% diagnostic_variables(jn))170 DO jn=1,size(model%interior_diagnostic_variables) 171 CALL write_variable_xml(xml_unit,model%interior_diagnostic_variables(jn),3) 172 CALL write_25hourm_xml(xml_unit,model%interior_diagnostic_variables(jn),3) 173 CALL write_tmb_xml(xml_unit,model%interior_diagnostic_variables(jn)) 161 174 END DO 162 175 DO jn=1,size(model%horizontal_diagnostic_variables) 163 176 CALL write_variable_xml(xml_unit,model%horizontal_diagnostic_variables(jn)) 164 177 CALL write_25hourm_xml(xml_unit,model%horizontal_diagnostic_variables(jn)) 178 END DO 179 DO jn=1,size(model%interior_state_variables) 180 WRITE (xml_unit,'(A)') ' <field id="'//TRIM(model%interior_state_variables(jn)%name)// & 181 & '_VINT" long_name="depth-integrated '//TRIM(model%interior_state_variables(jn)%long_name)// & 182 & '" unit="'//TRIM(model%interior_state_variables(jn)%units)//'*m" default_value="0.0" />' 183 END DO 184 DO jn=1,size(model%interior_diagnostic_variables) 185 WRITE (xml_unit,'(A)') ' <field id="'//TRIM(model%interior_diagnostic_variables(jn)%name)// & 186 & '_VINT" long_name="depth-integrated '//TRIM(model%interior_diagnostic_variables(jn)%long_name)// & 187 & '" unit="'//TRIM(model%interior_diagnostic_variables(jn)%units)//'*m" default_value="0.0" />' 165 188 END DO 166 189 WRITE (xml_unit,1000) ' </field_group>' … … 195 218 1000 FORMAT (A) 196 219 197 END SUBROUTINE nemo_fabm_ init220 END SUBROUTINE nemo_fabm_configure 198 221 199 222 SUBROUTINE write_variable_xml(xml_unit,variable,flag_grid_ref) 200 223 INTEGER,INTENT(IN) :: xml_unit 201 224 INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 202 CLASS (type_ external_variable),INTENT(IN) :: variable225 CLASS (type_fabm_variable),INTENT(IN) :: variable 203 226 204 227 CHARACTER(LEN=20) :: missing_value,string_dimensions … … 233 256 INTEGER,INTENT(IN) :: xml_unit 234 257 INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 235 CLASS (type_ external_variable),INTENT(IN) :: variable258 CLASS (type_fabm_variable),INTENT(IN) :: variable 236 259 237 260 CHARACTER(LEN=20) :: missing_value,string_dimensions … … 265 288 SUBROUTINE write_tmb_xml(xml_unit,variable) 266 289 INTEGER,INTENT(IN) :: xml_unit 267 CLASS (type_ external_variable),INTENT(IN) :: variable290 CLASS (type_fabm_variable),INTENT(IN) :: variable 268 291 269 292 CHARACTER(LEN=20) :: missing_value … … 279 302 INTEGER,INTENT(IN) :: xml_unit 280 303 INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 281 CLASS (type_ external_variable),INTENT(IN) :: variable304 CLASS (type_fabm_variable),INTENT(IN) :: variable 282 305 283 306 INTEGER :: number_dimensions,i … … 383 406 !! ** Purpose : initialization for FABM model 384 407 !! 385 !! ** Method : - Read the namcfc namelist and check the parameter values408 !! ** Method : - Allocate FABM arrays, configure domain, send data 386 409 !!---------------------------------------------------------------------- 387 410 #if defined key_git_version 388 TYPE (type_version), POINTER :: version411 TYPE (type_version), POINTER :: version 389 412 #endif 390 413 INTEGER :: jn 391 392 ! ! Allocate FABM arrays393 IF( trc_sms_fabm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trc_ini_fabm: unable to allocate FABM arrays' )394 414 395 415 IF(lwp) WRITE(numout,*) … … 399 419 IF(lwp) WRITE(numout,*) ' NEMO version: ',git_commit_id,' (',git_branch_name,' branch)' 400 420 IF(lwp) WRITE(numout,*) ' FABM version: ',fabm_commit_id,' (',fabm_branch_name,' branch)' 401 #endif402 403 call fabm_initialize_library()404 #if defined key_git_version405 421 version => first_module_version 406 407 422 do while (associated(version)) 408 423 IF(lwp) WRITE(numout,*) ' '//trim(version%module_name)//' version: ',trim(version%version_string) … … 411 426 #endif 412 427 428 ! Allocate FABM arrays 429 IF(trc_sms_fabm_alloc() /= 0) CALL ctl_stop( 'STOP', 'trc_ini_fabm: unable to allocate FABM arrays' ) 430 413 431 ! Log mapping of FABM states: 414 432 IF (lwp) THEN 415 IF (jp_fabm .gt.0) WRITE(numout,*) " FABM tracers:"433 IF (jp_fabm > 0) WRITE(numout,*) " FABM tracers:" 416 434 DO jn=1,jp_fabm 417 WRITE(numout,*) " State",jn,":",trim(model% state_variables(jn)%name), &418 " (",trim(model% state_variables(jn)%long_name), &419 ") [",trim(model% state_variables(jn)%units),"]"420 END DO421 IF (jp_fabm_surface .gt.0) WRITE(numout,*) "FABM seasurface states:"435 WRITE(numout,*) " State",jn,":",trim(model%interior_state_variables(jn)%name), & 436 " (",trim(model%interior_state_variables(jn)%long_name), & 437 ") [",trim(model%interior_state_variables(jn)%units),"]" 438 END DO 439 IF (jp_fabm_surface > 0) WRITE(numout,*) "FABM seasurface states:" 422 440 DO jn=1,jp_fabm_surface 423 441 WRITE(numout,*) " State",jn,":",trim(model%surface_state_variables(jn)%name), & 424 442 " (",trim(model%surface_state_variables(jn)%long_name), & 425 443 ") [",trim(model%surface_state_variables(jn)%units),"]" 426 END DO427 IF (jp_fabm_bottom .gt.0) WRITE(numout,*) "FABM seafloor states:"444 END DO 445 IF (jp_fabm_bottom > 0) WRITE(numout,*) "FABM seafloor states:" 428 446 DO jn=1,jp_fabm_bottom 429 447 WRITE(numout,*) " State",jn,":",trim(model%bottom_state_variables(jn)%name), & 430 448 " (",trim(model%bottom_state_variables(jn)%long_name), & 431 449 ") [",trim(model%bottom_state_variables(jn)%units),"]" 432 END DO433 END IF450 END DO 451 END IF 434 452 435 453 ! Initialise variables required for Hemmings et al. (2008) … … 442 460 END SUBROUTINE trc_ini_fabm 443 461 462 SUBROUTINE nemo_fabm_driver_fatal_error(self, location, message) 463 CLASS (type_nemo_fabm_driver), INTENT(INOUT) :: self 464 CHARACTER(len=*), INTENT(IN) :: location, message 465 466 CALL ctl_stop('STOP', TRIM(location)//': '//TRIM(message)) 467 STOP 468 END SUBROUTINE 469 470 SUBROUTINE nemo_fabm_driver_log_message(self, message) 471 CLASS (type_nemo_fabm_driver), INTENT(INOUT) :: self 472 CHARACTER(len=*), INTENT(IN) :: message 473 474 IF(lwp) WRITE (numout,*) TRIM(message) 475 END SUBROUTINE 476 444 477 INTEGER FUNCTION fabm_state_index( state_name ) 445 478 !!---------------------------------------------------------------------- … … 453 486 IMPLICIT NONE 454 487 455 CHARACTER(LEN= 256), INTENT(IN) :: state_name456 457 INTEGER 488 CHARACTER(LEN=*), INTENT(IN) :: state_name 489 490 INTEGER :: jn 458 491 459 492 !!---------------------------------------------------------------------- … … 461 494 fabm_state_index = -1 462 495 DO jn=1,jp_fabm 463 IF (TRIM(model% state_variables(jn)%name) == TRIM(state_name)) THEN496 IF (TRIM(model%interior_state_variables(jn)%name) == TRIM(state_name)) THEN 464 497 fabm_state_index = jn 465 498 EXIT … … 485 518 IMPLICIT NONE 486 519 487 CHARACTER(LEN= 256), INTENT(IN) :: diag_name488 489 INTEGER 520 CHARACTER(LEN=*), INTENT(IN) :: diag_name 521 522 INTEGER :: jn 490 523 491 524 !!---------------------------------------------------------------------- 492 525 493 526 fabm_diag_index = -1 494 DO jn = 1, SIZE(model% diagnostic_variables)495 IF (TRIM(model% diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN527 DO jn = 1, SIZE(model%interior_diagnostic_variables) 528 IF (TRIM(model%interior_diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN 496 529 fabm_diag_index = jn 497 530 EXIT -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90
r10270 r13576 4 4 !! TOP : initialisation of some run parameters for FABM bio-model 5 5 !!====================================================================== 6 !! History : 2.0 ! 2007-12 (C. Ethe, G. Madec) Original code 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 … … 18 19 USE par_fabm 19 20 USE trcsms_fabm 20 21 21 22 22 IMPLICIT NONE … … 35 35 36 36 SUBROUTINE trc_nam_fabm 37 LOGICAL :: l_ext 38 INTEGER :: nmlunit, ios 39 NAMELIST/namfabm/ nn_adv 40 41 ! Read NEMO-FABM coupler settings from namfabm 42 nn_adv = 3 43 INQUIRE( FILE='namelist_fabm_ref', EXIST=l_ext ) 44 IF (l_ext) then 45 CALL ctl_opn( nmlunit, 'namelist_fabm_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE.) 46 READ(nmlunit, nml=namfabm, iostat=ios) 47 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namfabm in namelist_fabm_ref', .TRUE. ) 48 END IF 49 INQUIRE( FILE='namelist_fabm_cfg', EXIST=l_ext ) 50 IF (l_ext) then 51 CALL ctl_opn( nmlunit, 'namelist_fabm_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE.) 52 READ(nmlunit, nml=namfabm, iostat=ios) 53 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namfabm in namelist_fabm_cfg', .TRUE. ) 54 END IF 55 IF (nn_adv /= 1 .and. nn_adv /= 3) CALL ctl_stop( 'STOP', 'trc_ini_fabm: nn_adv must be 1 or 3.' ) 37 56 END SUBROUTINE trc_nam_fabm 38 57 39 SUBROUTINE trc_nam_fabm_override 58 SUBROUTINE trc_nam_fabm_override(sn_tracer) 59 TYPE(PTRACER), DIMENSION(jpmaxtrc), INTENT(INOUT) :: sn_tracer 60 40 61 INTEGER :: jn 62 CHARACTER(LEN=3) :: index 41 63 42 64 DO jn=1,jp_fabm 43 ctrcnm(jp_fabm_m1+jn) = model%state_variables(jn)%name 44 ctrcln(jp_fabm_m1+jn) = model%state_variables(jn)%long_name 45 ctrcun(jp_fabm_m1+jn) = model%state_variables(jn)%units 46 ln_trc_ini(jp_fabm_m1+jn) = .FALSE. 65 IF (sn_tracer(jn)%clsname /= 'NONAME' .AND. sn_tracer(jn)%clsname /= model%interior_state_variables(jn)%name) THEN 66 WRITE (index,'(i0)') jn 67 CALL ctl_stop('Tracer name mismatch in namtrc: '//TRIM(sn_tracer(jn)%clsname)//' found at sn_tracer('//TRIM(index)//') where '//TRIM(model%interior_state_variables(jn)%name)//' was expected.') 68 END IF 69 sn_tracer(jn)%clsname = TRIM(model%interior_state_variables(jn)%name) 70 sn_tracer(jn)%cllname = TRIM(model%interior_state_variables(jn)%long_name) 71 sn_tracer(jn)%clunit = TRIM(model%interior_state_variables(jn)%units) 72 sn_tracer(jn)%llinit = .FALSE. 47 73 END DO 48 74 END SUBROUTINE trc_nam_fabm_override 49 75 50 76 #else 51 77 !!---------------------------------------------------------------------- -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90
r10156 r13576 4 4 !! Read and write additional restart fields used by FABM 5 5 !!====================================================================== 6 !! History : 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 -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90
r10728 r13576 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_pgrow)%save = .TRUE. 580 model%interior_diagnostic_variables(jp_fabm_ploss)%save = .TRUE. 581 IF( ln_qsr_spec ) THEN 582 model%interior_diagnostic_variables(jp_fabm_swr)%save = .TRUE. 583 ENDIF 584 IF ( jp_fabm_kd490 /= -1 ) THEN 585 model%interior_diagnostic_variables(jp_fabm_kd490)%save = .TRUE. 586 ENDIF 587 IF ( jp_fabm_xeps /= -1 ) THEN 588 model%interior_diagnostic_variables(jp_fabm_xeps)%save = .TRUE. 589 ENDIF 590 DO jn=1,size(model%horizontal_diagnostic_variables) 591 model%horizontal_diagnostic_variables(jn)%save = iom_use(model%horizontal_diagnostic_variables(jn)%name) & 592 .or. iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h') 593 END DO 594 595 ! Check whether FABM has all required data 596 ! [after this, the save attribute of diagnostic variables can no longer change!] 597 CALL model%start() 598 599 started = .TRUE. 600 601 END SUBROUTINE 574 602 575 603 #else -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90
r10270 r13576 4 4 !! fabm : Output of FABM tracers 5 5 !!====================================================================== 6 !! History : 1.0 ! 2009-05 (C. Ethe) Original code 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_top && key_fabm && defined key_iomput … … 18 19 USE par_fabm 19 20 USE st2d_fabm 20 USE fabm, only: fabm_get_bulk_diagnostic_data, fabm_get_horizontal_diagnostic_data21 21 22 22 IMPLICIT NONE … … 31 31 MODULE PROCEDURE wri_fabm,wri_fabm_fl 32 32 END INTERFACE trc_wri_fabm 33 34 33 35 34 PUBLIC trc_wri_fabm … … 58 57 ! depth integrated 59 58 ! for strict budgetting write this out at end of timestep as an average between 'now' and 'after' at kt 60 DO jn = 1, jp_fabm 161 IF(ln_trdtrc (j n))THEN62 trpool(:,:,:) = 0.5 * ( trn(:,:,:,jp_fabm 0+jn-1)*fse3t_a(:,:,:) + &59 DO jn = 1, jp_fabm 60 IF(ln_trdtrc (jp_fabm_m1+jn))THEN 61 trpool(:,:,:) = 0.5 * ( trn(:,:,:,jp_fabm_m1+jn)*fse3t_a(:,:,:) + & 63 62 tr_temp(:,:,:,jn)*fse3t(:,:,:) ) 64 cltra = TRIM( model% state_variables(jn)%name )//"_e3t" ! depth integrated output63 cltra = TRIM( model%interior_state_variables(jn)%name )//"_e3t" ! depth integrated output 65 64 IF( kt == nittrc000 ) write(6,*)'output pool ',cltra 66 65 CALL iom_put( cltra, trpool) … … 80 79 !!--------------------------------------------------------------------- 81 80 INTEGER, INTENT( in ) :: kt 82 INTEGER :: jn 81 INTEGER :: jn, jk 82 REAL(wp), DIMENSION(jpi,jpj) :: vint 83 83 84 84 #if defined key_tracer_budget … … 90 90 #endif 91 91 DO jn = 1, jp_fabm 92 CALL iom_put( model%state_variables(jn)%name, trn(:,:,:,jp_fabm0+jn-1) ) 92 ! Save 3D field 93 CALL iom_put(model%interior_state_variables(jn)%name, trn(:,:,:,jp_fabm_m1+jn)) 94 95 ! Save depth integral if selected for output in XIOS 96 IF (iom_use(TRIM(model%interior_state_variables(jn)%name)//'_VINT')) THEN 97 vint = 0._wp 98 DO jk = 1, jpkm1 99 vint = vint + trn(:,:,jk,jp_fabm_m1+jn) * fse3t(:,:,jk) * tmask(:,:,jk) 100 END DO 101 CALL iom_put(TRIM(model%interior_state_variables(jn)%name)//'_VINT', vint) 102 END IF 93 103 END DO 94 104 DO jn = 1, jp_fabm_surface … … 99 109 END DO 100 110 101 ! write 3D diagnostics in the file102 ! ---------------------------------------103 DO jn = 1, size(model%diagnostic_variables)104 IF (model%diagnostic_variables(jn)%save) &105 CALL iom_put( model%diagnostic_variables(jn)%name, fabm_get_bulk_diagnostic_data(model,jn))106 END DO107 108 ! write 2D diagnostics in the file109 ! ---------------------------------------110 DO jn = 1, size(model%horizontal_diagnostic_variables)111 IF (model%horizontal_diagnostic_variables(jn)%save) &112 CALL iom_put( model%horizontal_diagnostic_variables(jn)%name, fabm_get_horizontal_diagnostic_data(model,jn))113 END DO114 !115 111 CALL trc_sms_fabm_check_mass 116 112 … … 121 117 !! Dummy module : No passive tracer 122 118 !!---------------------------------------------------------------------- 119 INTERFACE trc_wri_fabm 120 MODULE PROCEDURE wri_fabm,wri_fabm_fl 121 END INTERFACE trc_wri_fabm 122 123 123 PUBLIC trc_wri_fabm 124 CONTAINS 125 SUBROUTINE trc_wri_fabm ! Empty routine 126 END SUBROUTINE trc_wri_fabm 124 125 CONTAINS 126 127 SUBROUTINE wri_fabm_fl(kt,fl) 128 INTEGER, INTENT( in ) :: fl 129 INTEGER, INTENT( in ) :: kt 130 END SUBROUTINE wri_fabm_fl 131 132 SUBROUTINE wri_fabm(kt) ! Empty routine 133 INTEGER, INTENT( in ) :: kt 134 END SUBROUTINE wri_fabm 135 127 136 #endif 128 137 -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90
r12352 r13576 15 15 USE trc 16 16 USE par_fabm 17 USE fabm18 17 USE dom_oce 19 18 #if defined key_trdtrc && defined key_iomput … … 25 24 26 25 # include "domzgr_substitute.h90" 26 # include "vectopt_loop_substitute.h90" 27 27 28 28 PRIVATE … … 32 32 ! Work arrays for vertical advection (residual movement/sinking/floating) 33 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: w_ct 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: w_if35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: zwgt_if36 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: flux_if37 34 #if defined key_trdtrc && defined key_iomput 38 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:,:) :: tr_vmv … … 41 38 CONTAINS 42 39 43 SUBROUTINE compute_vertical_movement( kt )40 SUBROUTINE compute_vertical_movement( kt, method ) 44 41 !!---------------------------------------------------------------------- 45 42 !! *** compute_vertical_movement *** 46 43 !! 47 !! ** Purpose : compute vertical movement of FABM tracers 44 !! ** Purpose : compute vertical movement of FABM tracers through the water 45 !! (sinking/floating/active movement) 48 46 !! 49 !! ** Method : Sets additional vertical velocity field and computes 50 !! resulting advection using a conservative 3rd upwind 51 !! scheme with QUICKEST TVD limiter, based on GOTM 52 !! module adv_center.F90 (www.gotm.net). Currently assuming 53 !! zero flux at sea surface and sea floor. 47 !! ** Method : Retrieves additional vertical velocity field and applies 48 !! advection scheme. 54 49 !!---------------------------------------------------------------------- 55 50 ! 56 INTEGER, INTENT(in) :: kt ! ocean time-step index 57 INTEGER :: ji,jj,jk,jn,k_floor,n_iter,n_count 58 INTEGER,PARAMETER :: n_itermax=100 59 REAL(wp) :: cmax_no,z2dt 60 REAL(wp),DIMENSION(jpk) :: tr_it,tr_u,tr_d,tr_c,tr_slope,c_no,flux_lim 61 REAL(wp),DIMENSION(jpk) :: phi_lim,x_fac 51 INTEGER, INTENT(in) :: kt ! ocean time-step index 52 INTEGER, INTENT(in) :: method ! advection method (1: 1st order upstream, 3: 3rd order TVD with QUICKEST limiter) 53 54 INTEGER :: ji,jj,jk,jn,k_floor 55 REAL(wp) :: zwgt_if(1:jpkm1-1), dc(1:jpkm1), w_if(1:jpkm1-1), z2dt, h(1:jpkm1) 62 56 #if defined key_trdtrc 63 57 CHARACTER (len=20) :: cltra … … 69 63 70 64 IF( neuler == 0 .AND. kt == nittrc000 ) THEN 71 65 z2dt = rdt ! set time step size (Euler) 72 66 ELSE 73 67 z2dt = 2._wp * rdt ! set time step size (Leapfrog) 74 68 ENDIF 69 75 70 ! Compute interior vertical velocities and include them in source array. 76 DO jj= 1,jpj! j-loop77 ! Get vertical velocities at layer centres (entire 1:jpi,1:jpk slice).78 DO jk=1,jpk 79 CALL fabm_get_vertical_movement(model,1,jpi,jj,jk,w_ct(:,jk,:))71 DO jj=2,jpjm1 ! j-loop 72 ! Get vertical velocities at layer centres (entire i-k slice). 73 DO jk=1,jpkm1 74 CALL model%get_vertical_movement(fs_2,fs_jpim1,jj,jk,w_ct(:,jk,:)) 80 75 END DO 81 82 DO ji=1,jpi ! i-loop 76 DO ji=fs_2,fs_jpim1 ! i-loop 83 77 ! Only process this horizontal point (ji,jj) if number of layers exceeds 1 84 IF (mbkt(ji,jj)>1) THEN ! Level check85 k_floor=mbkt(ji,jj)78 k_floor = mbkt(ji,jj) 79 IF (k_floor > 1) THEN 86 80 ! Linearly interpolate to velocities at the interfaces between layers 87 81 ! Note: 88 ! - interface k sits between cell centre k and k -1,89 ! - k [1,jpk ] increases downwards82 ! - interface k sits between cell centre k and k+1 (k=0 for surface) 83 ! - k [1,jpkm1] increases downwards 90 84 ! - upward velocity is positive, downward velocity is negative 91 zwgt_if(1,:)=0._wp ! surface 92 w_if(1,:)=0._wp ! surface 93 zwgt_if(2:k_floor,:)=spread(& 94 fse3t(ji,jj,2:k_floor)/ (fse3t(ji,jj,1:k_floor-1)+fse3t(ji,jj,2:k_floor))& 95 ,2,jp_fabm) 96 w_if(2:k_floor,:) = zwgt_if(2:k_floor,:)*w_ct(ji,1:k_floor-1,:)& 97 +(1._wp-zwgt_if(1:k_floor-1,:))*w_ct(ji,2:k_floor,:) 98 zwgt_if(k_floor+1:,:)=0._wp ! sea floor and below 99 w_if(k_floor+1:,:)=0._wp ! sea floor and below 85 h(1:k_floor) = fse3t(ji,jj,1:k_floor) 86 zwgt_if(1:k_floor-1) = h(2:k_floor) / (h(1:k_floor-1) + h(2:k_floor)) 100 87 101 88 ! Advect: 102 89 DO jn=1,jp_fabm ! State loop 103 ! get maximum Courant number: 104 c_no(2:k_floor)=abs(w_if(2:k_floor,jn))*z2dt/ & 105 ( 0.5_wp*(fse3t(ji,jj,2:k_floor) + & 106 fse3t(ji,jj,1:k_floor-1)) ) 107 cmax_no=MAXVAL(c_no(2:k_floor)) 108 109 ! number of iterations: 110 n_iter=min(n_itermax,int(cmax_no)+1) 111 IF (ln_ctl.AND.(n_iter .gt. 1)) THEN 112 WRITE(numout,*) 'vertical_movement_fabm():' 113 WRITE(numout,*) ' Maximum Courant number is ',cmax_no,'.' 114 WRITE(numout,*) ' ',n_iter,' iterations used for vertical advection.' 115 ENDIF 116 117 ! effective Courant number: 118 c_no=c_no/n_iter 119 120 tr_it(1:k_floor)=trb(ji,jj,1:k_floor,jp_fabm_m1+jn) 121 DO n_count=1,n_iter ! Iterative loop 122 !Compute slope ratio 123 IF (k_floor.gt.2) THEN !More than 2 vertical wet points 124 IF (k_floor.gt.3) THEN 125 WHERE (w_if(3:k_floor-1,jn).ge.0._wp) !upward movement 126 tr_u(3:k_floor-1)=tr_it(4:k_floor) 127 tr_c(3:k_floor-1)=tr_it(3:k_floor-1) 128 tr_d(3:k_floor-1)=tr_it(2:k_floor-2) 129 ELSEWHERE !downward movement 130 tr_u(3:k_floor-1)=tr_it(1:k_floor-3) 131 tr_c(3:k_floor-1)=tr_it(2:k_floor-2) 132 tr_d(3:k_floor-1)=tr_it(3:k_floor-1) 133 ENDWHERE 134 ENDIF 135 IF (w_if(2,jn).ge.0._wp) THEN 136 tr_u(2)=tr_it(3) 137 tr_c(2)=tr_it(2) 138 tr_d(2)=tr_it(1) 139 ELSE 140 tr_u(2)=tr_it(1) 141 tr_c(2)=tr_it(1) 142 tr_d(2)=tr_it(2) 143 ENDIF 144 IF (w_if(k_floor,jn).ge.0._wp) THEN 145 tr_u(k_floor)=tr_it(k_floor) 146 tr_c(k_floor)=tr_it(k_floor) 147 tr_d(k_floor)=tr_it(k_floor-1) 148 ELSE 149 tr_u(k_floor)=tr_it(k_floor-2) 150 tr_c(k_floor)=tr_it(k_floor-1) 151 tr_d(k_floor)=tr_it(k_floor) 152 ENDIF 153 ELSE !only 2 vertical wet points, i.e. only 1 interface 154 IF (w_if(k_floor,jn).ge.0._wp) THEN 155 tr_u(2)=tr_it(2) 156 tr_c(2)=tr_it(2) 157 tr_d(2)=tr_it(1) 158 ELSE 159 tr_u(2)=tr_it(1) 160 tr_c(2)=tr_it(1) 161 tr_d(2)=tr_it(2) 162 ENDIF 163 ENDIF 164 WHERE (abs(tr_d(2:k_floor)-tr_c(2:k_floor)).gt.1.e-10_wp) 165 tr_slope(2:k_floor)= & 166 (tr_c(2:k_floor)-tr_u(2:k_floor))/ & 167 (tr_d(2:k_floor)-tr_c(2:k_floor)) 168 ELSEWHERE 169 tr_slope(2:k_floor)=SIGN(1._wp,w_if(2:k_floor,jn))* & 170 (tr_c(2:k_floor)-tr_u(2:k_floor))*1.e10_wp 171 ENDWHERE 172 173 !QUICKEST flux limiter: 174 x_fac(2:k_floor)=(1._wp-2._wp*c_no(2:k_floor))/6._wp 175 phi_lim(2:k_floor)=(0.5_wp+x_fac(2:k_floor)) + & 176 (0.5_wp-x_Fac(2:k_floor))*tr_slope(2:k_floor) 177 flux_lim(2:k_floor)=max( 0._wp, & 178 min( phi_lim(2:k_floor),2._wp/(1._wp-c_no(2:k_floor)), & 179 2._wp*tr_slope(2:k_floor)/(c_no(2:k_floor)+1.e-10_wp)) ) 180 181 ! Compute limited flux: 182 flux_if(2:k_floor,jn) = w_if(2:k_floor,jn)* & 183 ( tr_c(2:k_floor) + & 184 0.5_wp*flux_lim(2:k_floor)*(1._wp-c_no(2:k_floor))* & 185 (tr_d(2:k_floor)-tr_c(2:k_floor)) ) 186 187 ! Compute pseudo update for trend aggregation: 188 tr_it(1:k_floor-1) = tr_it(1:k_floor-1) + & 189 z2dt/float(n_iter)/fse3t(ji,jj,1:k_floor-1)* & 190 flux_if(2:k_floor,jn) 191 tr_it(2:k_floor) = tr_it(2:k_floor) - & 192 z2dt/float(n_iter)/fse3t(ji,jj,2:k_floor)* & 193 flux_if(2:k_floor,jn) 194 195 ENDDO ! Iterative loop 196 197 ! Estimate rate of change from pseudo state updates (source 198 ! splitting): 199 tra(ji,jj,1:k_floor,jp_fabm_m1+jn) = & 200 tra(ji,jj,1:k_floor,jp_fabm_m1+jn) + & 201 (tr_it(1:k_floor) - trb(ji,jj,1:k_floor,jp_fabm_m1+jn))/z2dt 202 #if defined key_trdtrc && defined key_iomput 203 IF( lk_trdtrc .AND. ln_trdtrc( jp_fabm_m1+jn ) ) THEN 204 tr_vmv(ji,jj,1:k_floor,jn)=(tr_it(1:k_floor) - trb(ji,jj,1:k_floor,jn))/z2dt 90 IF (ALL(w_ct(ji,1:k_floor,jn) == 0._wp)) CYCLE 91 92 ! Compute velocities at interfaces 93 w_if(1:k_floor-1) = zwgt_if(1:k_floor-1) * w_ct(ji,1:k_floor-1,jn) + (1._wp - zwgt_if(1:k_floor-1)) * w_ct(ji,2:k_floor,jn) 94 95 ! Compute change (per volume) due to vertical movement per layer 96 IF (method == 1) THEN 97 CALL advect_1(k_floor, trn(ji,jj,1:k_floor,jp_fabm_m1+jn), w_if(1:k_floor-1), h(1:k_floor), z2dt, dc(1:k_floor)) 98 ELSE 99 CALL advect_3(k_floor, trb(ji,jj,1:k_floor,jp_fabm_m1+jn), w_if(1:k_floor-1), h(1:k_floor), z2dt, dc(1:k_floor)) 205 100 END IF 206 #endif 207 ENDDO ! State loop 101 102 ! Incorporate change due to vertical movement in sources-sinks 103 tra(ji,jj,1:k_floor,jp_fabm_m1+jn) = tra(ji,jj,1:k_floor,jp_fabm_m1+jn) + dc(1:k_floor) 104 105 #if defined key_trdtrc && defined key_iomput 106 ! Store change due to vertical movement as diagnostic 107 IF( lk_trdtrc .AND. ln_trdtrc( jp_fabm_m1+jn)) tr_vmv(ji,jj,1:k_floor,jn) = dc(1:k_floor) 108 #endif 109 END DO ! State loop 208 110 END IF ! Level check 209 111 END DO ! i-loop 210 112 END DO ! j-loop 113 211 114 #if defined key_trdtrc && defined key_iomput 212 115 DO jn=1,jp_fabm ! State loop … … 220 123 END SUBROUTINE compute_vertical_movement 221 124 125 SUBROUTINE advect_1(nk, c, w, h, dt, trend) 126 INTEGER, INTENT(IN) :: nk 127 REAL(wp), INTENT(IN) :: c(1:nk) 128 REAL(wp), INTENT(IN) :: w(1:nk-1) 129 REAL(wp), INTENT(IN) :: h(1:nk) 130 REAL(wp), INTENT(IN) :: dt 131 REAL(wp), INTENT(OUT) :: trend(1:nk) 132 133 REAL(wp) :: flux(0:nk) 134 INTEGER :: jk 135 ! Compute fluxes (per surface area) over at interfaces (remember: positive for upwards) 136 flux(0) = 0._wp 137 DO jk=1,nk-1 ! k-loop 138 IF (w(jk) > 0) THEN 139 ! Upward movement (source layer is jk+1) 140 flux(jk) = min(w(jk), h(jk+1)/dt) * c(jk+1) 141 ELSE 142 ! Downward movement (source layer is jk) 143 flux(jk) = max(w(jk), -h(jk)/dt) * c(jk) 144 END IF 145 END DO 146 flux(nk) = 0._wp 147 trend = (flux(1:nk) - flux(0:nk-1)) / h 148 END SUBROUTINE 149 150 SUBROUTINE advect_3(nk, c_old, w, h, dt, trend) 151 INTEGER, INTENT(IN) :: nk 152 REAL(wp), INTENT(IN) :: c_old(1:nk) 153 REAL(wp), INTENT(IN) :: w(1:nk-1) 154 REAL(wp), INTENT(IN) :: h(1:nk) 155 REAL(wp), INTENT(IN) :: dt 156 REAL(wp), INTENT(OUT) :: trend(1:nk) 157 158 INTEGER, PARAMETER :: n_itermax=100 159 REAL(wp) :: cmax_no 160 REAL(wp) :: cfl(1:nk-1) 161 INTEGER :: n_iter, n_count, jk 162 REAL(wp) :: c(1:nk) 163 REAL(wp) :: tr_u(1:nk-1) 164 REAL(wp) :: tr_c(1:nk-1) 165 REAL(wp) :: tr_d(1:nk-1) 166 REAL(wp) :: delta_tr_u(1:nk-1) 167 REAL(wp) :: delta_tr(1:nk-1) 168 REAL(wp) :: ratio(1:nk-1) 169 REAL(wp) :: x_fac(1:nk-1) 170 REAL(wp) :: phi_lim(1:nk-1) 171 REAL(wp) :: limiter(1:nk-1) 172 REAL(wp) :: flux_if(1:nk-1) 173 174 c(:) = c_old(:) 175 176 ! get maximum Courant number: 177 cfl = ABS(w) * dt / (0.5_wp * (h(2:nk) + h(1:nk-1))) 178 cmax_no = MAXVAL(cfl) 179 180 ! number of iterations: 181 n_iter = MIN(n_itermax, INT(cmax_no) + 1) 182 IF (ln_ctl.AND.(n_iter .gt. 1)) THEN 183 WRITE(numout,*) 'compute_vertical_movement::advect_3():' 184 WRITE(numout,*) ' Maximum Courant number is ',cmax_no,'.' 185 WRITE(numout,*) ' ',n_iter,' iterations used for vertical advection.' 186 ENDIF 187 188 ! effective Courant number: 189 cfl = cfl/n_iter 190 191 DO n_count=1,n_iter ! Iterative loop 192 ! Determine tracer concentration at 1.5 upstream (tr_u), 0.5 upstream (tr_c), 0.5 downstream (tr_d) from interface 193 IF (nk.gt.2) THEN 194 ! More than 2 vertical wet points 195 IF (nk.gt.3) THEN 196 WHERE (w(2:nk-2).ge.0._wp) 197 !upward movement 198 tr_u(2:nk-2)=c(4:nk) 199 tr_c(2:nk-2)=c(3:nk-1) 200 tr_d(2:nk-2)=c(2:nk-2) 201 ELSEWHERE 202 ! downward movement 203 tr_u(2:nk-2)=c(1:nk-3) 204 tr_c(2:nk-2)=c(2:nk-2) 205 tr_d(2:nk-2)=c(3:nk-1) 206 ENDWHERE 207 ENDIF 208 209 ! Interface between surface layer and the next 210 IF (w(1).ge.0._wp) THEN 211 ! upward movement 212 tr_u(1)=c(3) 213 tr_c(1)=c(2) 214 tr_d(1)=c(1) 215 ELSE 216 ! downward movement 217 tr_u(1)=c(1) 218 tr_c(1)=c(1) 219 tr_d(1)=c(2) 220 ENDIF 221 222 ! Interface between bottom layer and the previous 223 IF (w(nk-1).ge.0._wp) THEN 224 ! upward movement 225 tr_u(nk-1)=c(nk) 226 tr_c(nk-1)=c(nk) 227 tr_d(nk-1)=c(nk-1) 228 ELSE 229 ! downward movement 230 tr_u(nk-1)=c(nk-2) 231 tr_c(nk-1)=c(nk-1) 232 tr_d(nk-1)=c(nk) 233 ENDIF 234 ELSE 235 ! only 2 vertical wet points, i.e. only 1 interface 236 IF (w(1).ge.0._wp) THEN 237 ! upward movement 238 tr_u(1)=c(2) 239 tr_c(1)=c(2) 240 tr_d(1)=c(1) 241 ELSE 242 ! downward movement 243 tr_u(1)=c(1) 244 tr_c(1)=c(1) 245 tr_d(1)=c(2) 246 ENDIF 247 ENDIF 248 249 delta_tr_u = tr_c - tr_u 250 delta_tr = tr_d - tr_c 251 WHERE (delta_tr * delta_tr_u > 0._wp) 252 ! Monotonic function over tr_u, tr_c, r_d 253 254 ! Compute slope ratio 255 ratio = delta_tr_u / delta_tr 256 257 ! QUICKEST flux limiter 258 x_fac = (1._wp - 2._wp * cfl) / 6._wp 259 phi_lim = (0.5_wp + x_fac) + (0.5_wp - x_fac) * ratio 260 limiter = MIN(phi_lim, 2._wp / (1._wp - cfl), 2._wp * ratio / (cfl + 1.e-10_wp)) 261 262 ! Compute limited flux 263 flux_if = w * (tr_c + 0.5_wp * limiter * (1._wp - cfl) * delta_tr) 264 ELSEWHERE 265 ! Non-monotonic, use 1st order upstream 266 flux_if = w * tr_c 267 ENDWHERE 268 269 ! Compute pseudo update for trend aggregation: 270 c(1:nk-1) = c(1:nk-1) + dt / real(n_iter, wp) / h(1:nk-1) * flux_if 271 c(2:nk) = c(2:nk) - dt / real(n_iter, wp) / h(2:nk) * flux_if 272 273 ENDDO ! Iterative loop 274 275 ! Estimate rate of change from pseudo state updates (source splitting): 276 trend = (c - c_old) / dt 277 END SUBROUTINE 278 222 279 #endif 223 280 END MODULE -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/trcini.F90
r12432 r13576 76 76 ! +++>>> FABM 77 77 ! Allow FABM to update numbers of biogeochemical tracers, diagnostics (jptra etc.) 78 IF( lk_fabm ) CALL nemo_fabm_ init78 IF( lk_fabm ) CALL nemo_fabm_configure 79 79 ! FABM <<<+++ 80 80 … … 123 123 ! Initialisation of tracers Initial Conditions 124 124 IF( ln_trcdta ) CALL trc_dta_init(jptra) 125 126 125 127 126 IF( ln_rsttr ) THEN … … 162 161 ! FABM +++>>> 163 162 ! Initialisation of FABM diagnostics and tracer boundary conditions (so that you can use initial condition as boundary) 164 IF( lk_fabm ) THEN 165 wndm=0._wp !uninitiased field at this point 166 qsr=0._wp !uninitiased field at this point 167 CALL compute_fabm ! only needed to set-up diagnostics 168 CALL trc_bc_init(jptra) 169 ENDIF 163 IF( lk_fabm ) CALL trc_bc_init(jptra) 170 164 ! FABM <<<+++ 171 165 172 166 tra(:,:,:,:) = 0._wp 173 167 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav ) & ! Partial steps: before horizontal gradient of passive -
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/trcnam.F90
r10162 r13576 325 325 sn_tracer(:)%llcbc = .FALSE. 326 326 sn_tracer(:)%llcbc = .FALSE. 327 sn_tracer(:)%clsname = 'NONAME' 327 328 #endif 328 329 … … 335 336 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc in configuration namelist', lwp ) 336 337 IF(lwm) WRITE ( numont, namtrc ) 338 339 ! +++>>> FABM 340 if (lk_fabm) CALL trc_nam_fabm_override(sn_tracer) 341 ! FABM <<<+++ 337 342 338 343 DO jn = 1, jptra … … 354 359 END DO 355 360 356 ! +++>>> FABM357 if (lk_fabm) CALL trc_nam_fabm_override358 ! FABM <<<+++359 361 END SUBROUTINE trc_nam_trc 360 362 … … 380 382 INTEGER :: ios ! Local integer output status for namelist read 381 383 !!--------------------------------------------------------------------- 382 383 IF(lwp) WRITE(numout,*)384 IF(lwp) WRITE(numout,*) 'trc_nam_dia : read the passive tracer diagnostics options'385 IF(lwp) WRITE(numout,*) '~~~~~~~'386 384 387 385 IF(lwp) WRITE(numout,*)
Note: See TracChangeset
for help on using the changeset viewer.