Changeset 13241
- Timestamp:
- 2020-07-03T14:42:49+02:00 (4 years ago)
- Location:
- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO
- Files:
-
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
r10728 r13241 1082 1082 #elif defined key_fabm 1083 1083 totalk_tavg(:,:,:) = totalk_tavg(:,:,:) + & 1084 & fabm_get_interior_diagnostic_data(model,jp_fabm_o3ta) / pnumtimes_tavg1084 & model%get_interior_diagnostic_data(jp_fabm_o3ta) / pnumtimes_tavg 1085 1085 totalk_tavg(:,:,:) = totalk_tavg(:,:,:) * tmask(:,:,:) 1086 1086 #endif … … 1127 1127 cchl_p_tavg(:,:,:) = cchl_p(:,:,:) 1128 1128 #elif defined key_fabm 1129 totalk_tavg(:,:,:) = fabm_get_interior_diagnostic_data(model,jp_fabm_o3ta)1129 totalk_tavg(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_o3ta) 1130 1130 totalk_tavg(:,:,:) = totalk_tavg(:,:,:) * tmask(:,:,:) 1131 1131 #endif -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
r10390 r13241 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 fabm_3d_25h(:,:,:,jn) = model%get_interior_diagnostic_data(jn) 214 212 END DO 215 213 DO jn = 1, jp_fabm_surface … … 220 218 END DO 221 219 DO jn = 1, jp_fabm_2d 222 fabm_2d_25h(:,:,jn) = fabm_get_horizontal_diagnostic_data(model,jn)220 fabm_2d_25h(:,:,jn) = model%get_horizontal_diagnostic_data(jn) 223 221 END DO 224 222 #endif … … 327 325 END DO 328 326 DO jn = 1, jp_fabm_3d 329 fabm_3d_25h(:,:,:,jn) = fabm_3d_25h(:,:,:,jn) + fabm_get_interior_diagnostic_data(model,jn)327 fabm_3d_25h(:,:,:,jn) = fabm_3d_25h(:,:,:,jn) + model%get_interior_diagnostic_data(jn) 330 328 END DO 331 329 DO jn = 1, jp_fabm_surface … … 336 334 END DO 337 335 DO jn = 1, jp_fabm_2d 338 fabm_2d_25h(:,:,jn) = fabm_2d_25h(:,:,jn) + fabm_get_horizontal_diagnostic_data(model,jn)336 fabm_2d_25h(:,:,jn) = fabm_2d_25h(:,:,jn) + model%get_horizontal_diagnostic_data(jn) 339 337 END DO 340 338 #endif … … 401 399 DO jn = 1, jp_fabm 402 400 zw3d(:,:,:) = fabm_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 403 CALL iom_put( TRIM(model% state_variables(jn)%name)//"25h", zw3d )401 CALL iom_put( TRIM(model%interior_state_variables(jn)%name)//"25h", zw3d ) 404 402 END DO 405 403 DO jn = 1, jp_fabm_3d 406 404 zw3d(:,:,:) = fabm_3d_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 407 CALL iom_put( TRIM(model% diagnostic_variables(jn)%name)//"25h", zw3d )405 CALL iom_put( TRIM(model%interior_diagnostic_variables(jn)%name)//"25h", zw3d ) 408 406 END DO 409 407 DO jn = 1, jp_fabm_surface … … 468 466 END DO 469 467 DO jn = 1, jp_fabm_3d 470 fabm_3d_25h(:,:,:,jn) = fabm_get_interior_diagnostic_data(model,jn)468 fabm_3d_25h(:,:,:,jn) = model%get_interior_diagnostic_data(jn) 471 469 END DO 472 470 DO jn = 1, jp_fabm_surface … … 477 475 END DO 478 476 DO jn = 1, jp_fabm_2d 479 fabm_2d_25h(:,:,jn) = fabm_get_horizontal_diagnostic_data(model,jn)477 fabm_2d_25h(:,:,jn) = model%get_horizontal_diagnostic_data(jn) 480 478 END DO 481 479 #endif -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90
r10390 r13241 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 CALL dia_calctmb( model%get_interior_diagnostic_data(jn), zwtmb ) 179 CALL iom_put( "top_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,1) ) 180 CALL iom_put( "mid_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,2) ) 181 CALL iom_put( "bot_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,3) ) 183 182 END DO 184 183 #endif -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r8059 r13241 45 45 LOGICAL , PUBLIC :: ln_qsr_bio !: bio-model light absorption flag 46 46 LOGICAL , PUBLIC :: ln_qsr_ice !: light penetration for ice-model LIM3 (clem) 47 LOGICAL , PUBLIC :: ln_qsr_spec !: spectral model heating from ERSEM 47 48 INTEGER , PUBLIC :: nn_chldta !: use Chlorophyll data (=1) or not (=0) 48 49 INTEGER , PUBLIC :: nn_kd490dta !: use kd490dta data (=1) or not (=0) … … 150 151 151 152 ! ! ============================================== ! 152 IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! 153 IF( ln_qsr_spec ) THEN ! ERSEM spectral heating ! 154 ! ! ============================================== ! 155 DO jk = 1, jpkm1 156 qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) ) 157 END DO 158 ! Add to the general trend 159 DO jk = 1, jpkm1 160 DO jj = 2, jpjm1 161 DO ji = fs_2, fs_jpim1 ! vector opt. 162 z1_e3t = zfact / fse3t(ji,jj,jk) 163 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 164 END DO 165 END DO 166 END DO 167 CALL iom_put( 'qsr3d', etot3 ) ! Shortwave Radiation 3D distribution 168 IF ( ln_qsr_ice ) THEN 169 DO jj = 1, jpj 170 DO ji = 1, jpi 171 IF ( qsr(ji,jj) /= 0._wp ) THEN 172 fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) ) 173 ELSE 174 fraqsr_1lev(ji,jj) = 1. 175 ENDIF 176 END DO 177 END DO 178 ENDIF 179 ! 180 181 ! ! ============================================== ! 182 ELSEIF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! 153 183 ! ! ============================================== ! 154 184 DO jk = 1, jpkm1 … … 423 453 !! 424 454 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_kd490dta455 & ln_qsr_spec, nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta 426 456 !!---------------------------------------------------------------------- 427 457 … … 451 481 WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd 452 482 WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio 483 WRITE(numout,*) ' ERSEM spectral heating model ln_qsr_spec= ', ln_qsr_spec 453 484 WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice 454 485 WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta … … 470 501 IF( ln_qsr_2bd ) ioptio = ioptio + 1 471 502 IF( ln_qsr_bio ) ioptio = ioptio + 1 503 IF( ln_qsr_spec ) ioptio = ioptio + 1 472 504 IF( nn_kd490dta == 1 ) ioptio = ioptio + 1 473 505 ! … … 481 513 IF( ln_qsr_bio ) nqsr = 4 482 514 IF( nn_kd490dta == 1 ) nqsr = 5 515 IF( ln_qsr_spec ) nqsr = 6 483 516 ! 484 517 IF(lwp) THEN ! Print the choice … … 489 522 IF( nqsr == 4 ) WRITE(numout,*) ' bio-model light penetration' 490 523 IF( nqsr == 5 ) WRITE(numout,*) ' KD490 light penetration' 491 ENDIF 524 IF( nqsr == 6 ) WRITE(numout,*) ' ERSEM spectral light penetration' 525 ENDIF 526 #if ! defined key_fabm 527 ! 528 IF( nqsr == 6 ) THEN 529 CALL ctl_stop( 'ln_qsr_spec=.true. so trying to use ERSEM spectral light penetration', & 530 & 'but not running with ERSEM' ) 531 ENDIF 532 #endif 492 533 ! 493 534 ENDIF -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90
r10158 r13241 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_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90
r11545 r13241 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, & … … 37 41 jp_fabm_pgrow, jp_fabm_ploss 38 42 43 LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) :: lk_rad_fabm !: FABM negativity correction flag array 44 39 45 #if defined key_fabm 46 CLASS (type_fabm_model), POINTER :: model !FABM model instance 47 40 48 !!--------------------------------------------------------------------- 41 49 !! 'key_fabm' FABM tracers 42 50 !!--------------------------------------------------------------------- 43 51 LOGICAL, PUBLIC, PARAMETER :: lk_fabm = .TRUE. !: FABM flag 44 LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) :: lk_rad_fabm !: FABM negativity correction flag array45 52 #else 46 53 !!--------------------------------------------------------------------- -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90
r12352 r13241 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 … … 128 139 IF (lwp) THEN 129 140 ! write field_def_fabm.xml on lead process 130 OPEN(UNIT=xml_unit, FILE='field_def_fabm.xml',ACTION='WRITE',STATUS='REPLACE')141 OPEN(UNIT=xml_unit, FILE='field_def_fabm.xml', ACTION='WRITE', STATUS='REPLACE') 131 142 132 143 WRITE (xml_unit,1000) '<field_definition level="1" prec="4" operation="average" enabled=".TRUE." default_value="1.e20" >' … … 134 145 WRITE (xml_unit,1000) ' <field_group id="ptrc_T" grid_ref="grid_T_3D">' 135 146 DO jn=1,jp_fabm 136 CALL write_variable_xml(xml_unit,model% state_variables(jn))147 CALL write_variable_xml(xml_unit,model%interior_state_variables(jn)) 137 148 #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))149 CALL write_trends_xml(xml_unit,model%interior_state_variables(jn)) 150 #endif 151 CALL write_25hourm_xml(xml_unit,model%interior_state_variables(jn)) 152 CALL write_tmb_xml(xml_unit,model%interior_state_variables(jn)) 142 153 END DO 143 154 WRITE (xml_unit,1000) ' </field_group>' … … 155 166 156 167 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))168 DO jn=1,size(model%interior_diagnostic_variables) 169 CALL write_variable_xml(xml_unit,model%interior_diagnostic_variables(jn),3) 170 CALL write_25hourm_xml(xml_unit,model%interior_diagnostic_variables(jn),3) 171 CALL write_tmb_xml(xml_unit,model%interior_diagnostic_variables(jn)) 161 172 END DO 162 173 DO jn=1,size(model%horizontal_diagnostic_variables) 163 174 CALL write_variable_xml(xml_unit,model%horizontal_diagnostic_variables(jn)) 164 175 CALL write_25hourm_xml(xml_unit,model%horizontal_diagnostic_variables(jn)) 176 END DO 177 DO jn=1,size(model%interior_state_variables) 178 WRITE (xml_unit,'(A)') ' <field id="'//TRIM(model%interior_state_variables(jn)%name)// & 179 & '_VINT" long_name="depth-integrated '//TRIM(model%interior_state_variables(jn)%long_name)// & 180 & '" unit="'//TRIM(model%interior_state_variables(jn)%units)//'*m" default_value="0.0" />' 181 END DO 182 DO jn=1,size(model%interior_diagnostic_variables) 183 WRITE (xml_unit,'(A)') ' <field id="'//TRIM(model%interior_diagnostic_variables(jn)%name)// & 184 & '_VINT" long_name="depth-integrated '//TRIM(model%interior_diagnostic_variables(jn)%long_name)// & 185 & '" unit="'//TRIM(model%interior_diagnostic_variables(jn)%units)//'*m" default_value="0.0" />' 165 186 END DO 166 187 WRITE (xml_unit,1000) ' </field_group>' … … 195 216 1000 FORMAT (A) 196 217 197 END SUBROUTINE nemo_fabm_ init218 END SUBROUTINE nemo_fabm_configure 198 219 199 220 SUBROUTINE write_variable_xml(xml_unit,variable,flag_grid_ref) 200 221 INTEGER,INTENT(IN) :: xml_unit 201 222 INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 202 CLASS (type_ external_variable),INTENT(IN) :: variable223 CLASS (type_fabm_variable),INTENT(IN) :: variable 203 224 204 225 CHARACTER(LEN=20) :: missing_value,string_dimensions … … 233 254 INTEGER,INTENT(IN) :: xml_unit 234 255 INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 235 CLASS (type_ external_variable),INTENT(IN) :: variable256 CLASS (type_fabm_variable),INTENT(IN) :: variable 236 257 237 258 CHARACTER(LEN=20) :: missing_value,string_dimensions … … 265 286 SUBROUTINE write_tmb_xml(xml_unit,variable) 266 287 INTEGER,INTENT(IN) :: xml_unit 267 CLASS (type_ external_variable),INTENT(IN) :: variable288 CLASS (type_fabm_variable),INTENT(IN) :: variable 268 289 269 290 CHARACTER(LEN=20) :: missing_value … … 279 300 INTEGER,INTENT(IN) :: xml_unit 280 301 INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 281 CLASS (type_ external_variable),INTENT(IN) :: variable302 CLASS (type_fabm_variable),INTENT(IN) :: variable 282 303 283 304 INTEGER :: number_dimensions,i … … 383 404 !! ** Purpose : initialization for FABM model 384 405 !! 385 !! ** Method : - Read the namcfc namelist and check the parameter values406 !! ** Method : - Allocate FABM arrays, configure domain, send data 386 407 !!---------------------------------------------------------------------- 387 408 #if defined key_git_version 388 TYPE (type_version), POINTER :: version409 TYPE (type_version), POINTER :: version 389 410 #endif 390 411 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 412 395 413 IF(lwp) WRITE(numout,*) … … 399 417 IF(lwp) WRITE(numout,*) ' NEMO version: ',git_commit_id,' (',git_branch_name,' branch)' 400 418 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 419 version => first_module_version 406 407 420 do while (associated(version)) 408 421 IF(lwp) WRITE(numout,*) ' '//trim(version%module_name)//' version: ',trim(version%version_string) … … 411 424 #endif 412 425 426 ! Allocate FABM arrays 427 IF(trc_sms_fabm_alloc() /= 0) CALL ctl_stop( 'STOP', 'trc_ini_fabm: unable to allocate FABM arrays' ) 428 413 429 ! Log mapping of FABM states: 414 430 IF (lwp) THEN 415 IF (jp_fabm .gt.0) WRITE(numout,*) " FABM tracers:"431 IF (jp_fabm > 0) WRITE(numout,*) " FABM tracers:" 416 432 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:"433 WRITE(numout,*) " State",jn,":",trim(model%interior_state_variables(jn)%name), & 434 " (",trim(model%interior_state_variables(jn)%long_name), & 435 ") [",trim(model%interior_state_variables(jn)%units),"]" 436 END DO 437 IF (jp_fabm_surface > 0) WRITE(numout,*) "FABM seasurface states:" 422 438 DO jn=1,jp_fabm_surface 423 439 WRITE(numout,*) " State",jn,":",trim(model%surface_state_variables(jn)%name), & 424 440 " (",trim(model%surface_state_variables(jn)%long_name), & 425 441 ") [",trim(model%surface_state_variables(jn)%units),"]" 426 END DO427 IF (jp_fabm_bottom .gt.0) WRITE(numout,*) "FABM seafloor states:"442 END DO 443 IF (jp_fabm_bottom > 0) WRITE(numout,*) "FABM seafloor states:" 428 444 DO jn=1,jp_fabm_bottom 429 445 WRITE(numout,*) " State",jn,":",trim(model%bottom_state_variables(jn)%name), & 430 446 " (",trim(model%bottom_state_variables(jn)%long_name), & 431 447 ") [",trim(model%bottom_state_variables(jn)%units),"]" 432 END DO433 END IF448 END DO 449 END IF 434 450 435 451 ! Initialise variables required for Hemmings et al. (2008) … … 442 458 END SUBROUTINE trc_ini_fabm 443 459 460 SUBROUTINE nemo_fabm_driver_fatal_error(self, location, message) 461 CLASS (type_nemo_fabm_driver), INTENT(INOUT) :: self 462 CHARACTER(len=*), INTENT(IN) :: location, message 463 464 CALL ctl_stop('STOP', TRIM(location)//': '//TRIM(message)) 465 STOP 466 END SUBROUTINE 467 468 SUBROUTINE nemo_fabm_driver_log_message(self, message) 469 CLASS (type_nemo_fabm_driver), INTENT(INOUT) :: self 470 CHARACTER(len=*), INTENT(IN) :: message 471 472 IF(lwp) WRITE (numout,*) TRIM(message) 473 END SUBROUTINE 474 444 475 INTEGER FUNCTION fabm_state_index( state_name ) 445 476 !!---------------------------------------------------------------------- … … 461 492 fabm_state_index = -1 462 493 DO jn=1,jp_fabm 463 IF (TRIM(model% state_variables(jn)%name) == TRIM(state_name)) THEN494 IF (TRIM(model%interior_state_variables(jn)%name) == TRIM(state_name)) THEN 464 495 fabm_state_index = jn 465 496 EXIT … … 492 523 493 524 fabm_diag_index = -1 494 DO jn = 1, SIZE(model% diagnostic_variables)495 IF (TRIM(model% diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN525 DO jn = 1, SIZE(model%interior_diagnostic_variables) 526 IF (TRIM(model%interior_diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN 496 527 fabm_diag_index = jn 497 528 EXIT -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90
r10270 r13241 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_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90
r10156 r13241 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_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 -
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90
r10270 r13241 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_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90
r12352 r13241 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_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/trcini.F90
r12432 r13241 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_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/trcnam.F90
r10162 r13241 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.