New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 13241 – NEMO

Changeset 13241


Ignore:
Timestamp:
2020-07-03T14:42:49+02:00 (4 years ago)
Author:
dford
Message:

Update NEMO-FABM coupler for compatability with FABM v1.0.

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  
    10821082#elif defined key_fabm 
    10831083      totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) + & 
    1084          &                  fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) / pnumtimes_tavg 
     1084         &                  model%get_interior_diagnostic_data(jp_fabm_o3ta) / pnumtimes_tavg 
    10851085      totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:) 
    10861086#endif 
     
    11271127         cchl_p_tavg(:,:,:)     = cchl_p(:,:,:) 
    11281128#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) 
    11301130         totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:) 
    11311131#endif 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90

    r10390 r13241  
    2525   USE par_fabm 
    2626   USE st2d_fabm, ONLY: fabm_st2dn 
    27    USE fabm, ONLY: fabm_get_interior_diagnostic_data, & 
    28       &            fabm_get_horizontal_diagnostic_data 
    2927#endif 
    3028 
     
    211209      END DO 
    212210      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) 
    214212      END DO 
    215213      DO jn = 1, jp_fabm_surface 
     
    220218      END DO 
    221219      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) 
    223221      END DO 
    224222#endif 
     
    327325      END DO 
    328326      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) 
    330328      END DO 
    331329      DO jn = 1, jp_fabm_surface 
     
    336334      END DO 
    337335      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) 
    339337      END DO 
    340338#endif 
     
    401399            DO jn = 1, jp_fabm 
    402400               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  ) 
    404402            END DO 
    405403            DO jn = 1, jp_fabm_3d 
    406404               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  ) 
    408406            END DO 
    409407            DO jn = 1, jp_fabm_surface 
     
    468466      END DO 
    469467      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) 
    471469      END DO 
    472470      DO jn = 1, jp_fabm_surface 
     
    477475      END DO 
    478476      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) 
    480478      END DO 
    481479#endif 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90

    r10390 r13241  
    1414   USE trc, ONLY: trn 
    1515   USE par_fabm 
    16    USE fabm, ONLY: fabm_get_interior_diagnostic_data 
    1716#endif 
    1817 
     
    172171         DO jn = 1, jp_fabm 
    173172            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) ) 
    177176         END DO 
    178177         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) ) 
    183182         END DO 
    184183#endif 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r8059 r13241  
    4545   LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag 
    4646   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem) 
     47   LOGICAL , PUBLIC ::   ln_qsr_spec  !: spectral model heating from ERSEM 
    4748   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0) 
    4849   INTEGER , PUBLIC ::   nn_kd490dta  !: use kd490dta data (=1) or not (=0) 
     
    150151       
    151152      !                                           ! ============================================== ! 
    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  ! 
    153183         !                                        ! ============================================== ! 
    154184         DO jk = 1, jpkm1 
     
    423453      !! 
    424454      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_kd490dta 
     455         &                  ln_qsr_spec, nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta 
    426456      !!---------------------------------------------------------------------- 
    427457 
     
    451481         WRITE(numout,*) '      2 band               light penetration   ln_qsr_2bd = ', ln_qsr_2bd 
    452482         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 
    453484         WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice 
    454485         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)    nn_chldta  = ', nn_chldta 
     
    470501         IF( ln_qsr_2bd  )   ioptio = ioptio + 1 
    471502         IF( ln_qsr_bio  )   ioptio = ioptio + 1 
     503         IF( ln_qsr_spec )   ioptio = ioptio + 1 
    472504         IF( nn_kd490dta == 1 )   ioptio = ioptio + 1 
    473505         ! 
     
    481513         IF( ln_qsr_bio                      )   nqsr =  4 
    482514         IF( nn_kd490dta == 1                )   nqsr =  5 
     515         IF( ln_qsr_spec                     )   nqsr =  6 
    483516         ! 
    484517         IF(lwp) THEN                   ! Print the choice 
     
    489522            IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration' 
    490523            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 
    492533         ! 
    493534      ENDIF 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90

    r10158 r13241  
    1919   USE fldread 
    2020   USE par_fabm 
    21    USE fabm 
     21   USE fabm, only: type_fabm_horizontal_variable_id 
    2222 
    2323   IMPLICIT NONE 
     24 
     25#  include "vectopt_loop_substitute.h90" 
    2426 
    2527   PRIVATE 
     
    4042 
    4143   TYPE, PUBLIC, EXTENDS(type_input_variable) :: type_input_data 
    42       TYPE(type_horizontal_variable_id)  :: horizontal_id 
    43       TYPE(type_input_data), POINTER   :: next => null() 
     44      TYPE(type_fabm_horizontal_variable_id) :: horizontal_id 
     45      TYPE(type_input_data), POINTER         :: next => null() 
    4446   END TYPE 
    4547   TYPE (type_input_data), POINTER, PUBLIC :: first_input_data => NULL() 
     
    8789           ALLOCATE(input_data, STAT=ierr) 
    8890           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)) THEN 
     91           input_data%horizontal_id = model%get_horizontal_variable_id(name) 
     92           IF (.NOT.model%is_variable_used(input_data%horizontal_id)) THEN 
    9193              ! This variable was not found among FABM's horizontal variables (at least, those that are read by one or more FABM modules) 
    9294              CALL ctl_stop('STOP', 'inputs_fabm:initialize_inputs: variable "'//TRIM(name)//'" was not found among horizontal FABM variables.') 
     
    130132           ! within tracer field 
    131133           DO jn=1,jp_fabm 
    132              IF (TRIM(name) == TRIM(model%state_variables(jn)%name)) THEN 
     134             IF (TRIM(name) == TRIM(model%interior_state_variables(jn)%name)) THEN 
    133135               river_data%jp_pos = jp_fabm_m1+jn 
    134136             END IF 
     
    173175         ! Provide FABM with pointer to field that will receive prescribed data. 
    174176         ! 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) 
    176178         input_data => input_data%next 
    177179      END DO 
     
    226228#endif 
    227229        IF( kt == nit000 .OR. ( kt /= nit000 ) ) THEN 
    228             DO jj = 1, jpj 
    229               DO ji = 1, jpi 
     230            DO jj = 2, jpjm1 
     231              DO ji = fs_2, fs_jpim1 
    230232                ! convert units and divide by surface area 
    231233                ! 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  
    11MODULE par_fabm 
    22 
     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 
    38   USE fabm 
     9#endif 
    410 
    511   IMPLICIT NONE 
    6  
    7    TYPE (type_model) :: model !FABM model instance 
    812 
    913   INTEGER, PUBLIC :: jp_fabm0, jp_fabm1, jp_fabm, & 
     
    3741                      jp_fabm_pgrow, jp_fabm_ploss 
    3842 
     43   LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) ::   lk_rad_fabm !: FABM negativity correction flag array 
     44 
    3945#if defined key_fabm 
     46   CLASS (type_fabm_model), POINTER :: model !FABM model instance 
     47 
    4048   !!--------------------------------------------------------------------- 
    4149   !!   'key_fabm'                     FABM tracers 
    4250   !!--------------------------------------------------------------------- 
    4351   LOGICAL, PUBLIC, PARAMETER ::   lk_fabm     = .TRUE.   !: FABM flag  
    44    LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) ::   lk_rad_fabm !: FABM negativity correction flag array  
    4552#else 
    4653   !!--------------------------------------------------------------------- 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90

    r12352 r13241  
    44   !! TOP :   initialisation of the FABM tracers 
    55   !!====================================================================== 
    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 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_fabm 
     
    1718   USE par_fabm 
    1819   USE trcsms_fabm 
    19    USE fabm_config,ONLY: fabm_create_model_from_yaml_file 
    20    USE fabm,ONLY: type_external_variable, fabm_initialize_library 
    21    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, & 
    2223     type_input_variable,type_input_data,type_river_data, & 
    2324     first_input_data,first_river_data 
     
    3334 
    3435#if defined key_git_version 
    35 !#include "gitversion.h90" 
     36include "gitversion.h90" 
    3637   CHARACTER(len=*),parameter :: git_commit_id = _NEMO_COMMIT_ID_ 
    3738   CHARACTER(len=*),parameter :: git_branch_name = _NEMO_BRANCH_ 
     
    3940 
    4041   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 
    4249 
    4350   !!---------------------------------------------------------------------- 
     
    4855CONTAINS 
    4956 
    50    SUBROUTINE nemo_fabm_init() 
     57   SUBROUTINE nemo_fabm_configure() 
    5158      INTEGER :: jn 
    5259      INTEGER, PARAMETER :: xml_unit = 1979 
     
    5562      CLASS (type_input_variable),POINTER :: input_pointer 
    5663 
     64      ALLOCATE(type_nemo_fabm_driver::driver) 
     65 
    5766      ! 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) 
    6170      jp_fabm_bottom = size(model%bottom_state_variables) 
    6271      jp_fabm_surface = size(model%surface_state_variables) 
     
    6675      jptra = jptra + jp_fabm 
    6776      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) 
    6978      jpdia2d = jpdia2d + jp_fabm_2d 
    7079      jpdia3d = jpdia3d + jp_fabm_3d 
    7180      jpdiabio = jpdiabio + jp_fabm 
    7281 
    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. 
    7485      call initialize_inputs 
    7586 
     
    128139      IF (lwp) THEN 
    129140         ! 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') 
    131142 
    132143         WRITE (xml_unit,1000) '<field_definition level="1" prec="4" operation="average" enabled=".TRUE." default_value="1.e20" >' 
     
    134145         WRITE (xml_unit,1000) ' <field_group id="ptrc_T" grid_ref="grid_T_3D">' 
    135146         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)) 
    137148#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)) 
    142153         END DO 
    143154         WRITE (xml_unit,1000) ' </field_group>' 
     
    155166 
    156167         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)) 
    161172         END DO 
    162173         DO jn=1,size(model%horizontal_diagnostic_variables) 
    163174            CALL write_variable_xml(xml_unit,model%horizontal_diagnostic_variables(jn)) 
    164175            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" />' 
    165186         END DO 
    166187         WRITE (xml_unit,1000) ' </field_group>' 
     
    1952161000 FORMAT (A) 
    196217 
    197    END SUBROUTINE nemo_fabm_init 
     218   END SUBROUTINE nemo_fabm_configure 
    198219 
    199220   SUBROUTINE write_variable_xml(xml_unit,variable,flag_grid_ref) 
    200221      INTEGER,INTENT(IN) :: xml_unit 
    201222      INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 
    202       CLASS (type_external_variable),INTENT(IN) :: variable 
     223      CLASS (type_fabm_variable),INTENT(IN) :: variable 
    203224 
    204225      CHARACTER(LEN=20) :: missing_value,string_dimensions 
     
    233254      INTEGER,INTENT(IN) :: xml_unit 
    234255      INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 
    235       CLASS (type_external_variable),INTENT(IN) :: variable 
     256      CLASS (type_fabm_variable),INTENT(IN) :: variable 
    236257 
    237258      CHARACTER(LEN=20) :: missing_value,string_dimensions 
     
    265286   SUBROUTINE write_tmb_xml(xml_unit,variable) 
    266287      INTEGER,INTENT(IN) :: xml_unit 
    267       CLASS (type_external_variable),INTENT(IN) :: variable 
     288      CLASS (type_fabm_variable),INTENT(IN) :: variable 
    268289 
    269290      CHARACTER(LEN=20) :: missing_value 
     
    279300      INTEGER,INTENT(IN) :: xml_unit 
    280301      INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 
    281       CLASS (type_external_variable),INTENT(IN) :: variable 
     302      CLASS (type_fabm_variable),INTENT(IN) :: variable 
    282303 
    283304      INTEGER :: number_dimensions,i 
     
    383404      !! ** Purpose :   initialization for FABM model 
    384405      !! 
    385       !! ** Method  : - Read the namcfc namelist and check the parameter values 
     406      !! ** Method  : - Allocate FABM arrays, configure domain, send data 
    386407      !!---------------------------------------------------------------------- 
    387408#if defined key_git_version 
    388       TYPE (type_version),POINTER :: version 
     409      TYPE (type_version), POINTER :: version 
    389410#endif 
    390411      INTEGER :: jn 
    391  
    392       !                       ! Allocate FABM arrays 
    393       IF( trc_sms_fabm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trc_ini_fabm: unable to allocate FABM arrays' ) 
    394412 
    395413      IF(lwp) WRITE(numout,*) 
     
    399417      IF(lwp) WRITE(numout,*) ' NEMO version:   ',git_commit_id,' (',git_branch_name,' branch)' 
    400418      IF(lwp) WRITE(numout,*) ' FABM version:   ',fabm_commit_id,' (',fabm_branch_name,' branch)' 
    401 #endif 
    402  
    403       call fabm_initialize_library() 
    404 #if defined key_git_version 
    405419      version => first_module_version 
    406  
    407420      do while (associated(version)) 
    408421         IF(lwp) WRITE(numout,*)  ' '//trim(version%module_name)//' version:   ',trim(version%version_string) 
     
    411424#endif 
    412425 
     426      ! Allocate FABM arrays 
     427      IF(trc_sms_fabm_alloc() /= 0) CALL ctl_stop( 'STOP', 'trc_ini_fabm: unable to allocate FABM arrays' ) 
     428 
    413429      ! Log mapping of FABM states: 
    414430      IF (lwp) THEN 
    415          IF (jp_fabm.gt.0) WRITE(numout,*) " FABM tracers:" 
     431         IF (jp_fabm > 0) WRITE(numout,*) " FABM tracers:" 
    416432         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          ENDDO 
    421          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:" 
    422438         DO jn=1,jp_fabm_surface 
    423439            WRITE(numout,*) "   State",jn,":",trim(model%surface_state_variables(jn)%name), & 
    424440               " (",trim(model%surface_state_variables(jn)%long_name), & 
    425441               ") [",trim(model%surface_state_variables(jn)%units),"]" 
    426          ENDDO 
    427          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:" 
    428444         DO jn=1,jp_fabm_bottom 
    429445            WRITE(numout,*) "   State",jn,":",trim(model%bottom_state_variables(jn)%name), & 
    430446               " (",trim(model%bottom_state_variables(jn)%long_name), & 
    431447               ") [",trim(model%bottom_state_variables(jn)%units),"]" 
    432          ENDDO 
    433       ENDIF 
     448         END DO 
     449      END IF 
    434450       
    435451      ! Initialise variables required for Hemmings et al. (2008) 
     
    442458   END SUBROUTINE trc_ini_fabm 
    443459 
     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 
    444475   INTEGER FUNCTION fabm_state_index( state_name ) 
    445476      !!---------------------------------------------------------------------- 
     
    461492      fabm_state_index = -1 
    462493      DO jn=1,jp_fabm 
    463          IF (TRIM(model%state_variables(jn)%name) == TRIM(state_name)) THEN 
     494         IF (TRIM(model%interior_state_variables(jn)%name) == TRIM(state_name)) THEN 
    464495            fabm_state_index = jn 
    465496            EXIT 
     
    492523       
    493524      fabm_diag_index = -1 
    494       DO jn = 1, SIZE(model%diagnostic_variables) 
    495          IF (TRIM(model%diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN 
     525      DO jn = 1, SIZE(model%interior_diagnostic_variables) 
     526         IF (TRIM(model%interior_diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN 
    496527            fabm_diag_index = jn 
    497528            EXIT 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90

    r10270 r13241  
    44   !! TOP :   initialisation of some run parameters for FABM bio-model 
    55   !!====================================================================== 
    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 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_fabm 
     
    1819   USE par_fabm 
    1920   USE trcsms_fabm 
    20  
    2121 
    2222   IMPLICIT NONE 
     
    3535 
    3636   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.' ) 
    3756   END SUBROUTINE trc_nam_fabm 
    3857 
    39    SUBROUTINE trc_nam_fabm_override 
     58   SUBROUTINE trc_nam_fabm_override(sn_tracer) 
     59      TYPE(PTRACER), DIMENSION(jpmaxtrc), INTENT(INOUT) :: sn_tracer 
     60 
    4061      INTEGER :: jn 
     62      CHARACTER(LEN=3) :: index 
    4163 
    4264      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. 
    4773      END DO 
    4874   END SUBROUTINE trc_nam_fabm_override 
    49     
     75 
    5076#else 
    5177   !!---------------------------------------------------------------------- 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90

    r10156 r13241  
    44   !! Read and write additional restart fields used by FABM 
    55   !!====================================================================== 
    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 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_fabm 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90

    r10728 r13241  
    55   !!====================================================================== 
    66   !! History :   1.0  !  2015-04  (PML) Original code 
     7   !! History :   1.1  !  2020-06  (PML) Update to FABM 1.0, improved performance 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_fabm 
     
    1920   USE trd_oce 
    2021   USE trdtrc 
     22   USE traqsr,only: ln_qsr_spec 
    2123#if defined key_trdtrc && defined key_iomput 
    2224   USE trdtrc_oce 
     
    2426 
    2527   USE oce, only: tsn  ! Needed? 
    26    USE sbc_oce, only: lk_oasis 
     28   USE sbc_oce, only: lk_oasis,fr_i 
    2729   USE dom_oce 
    2830   USE zdf_oce 
    29    !USE iom 
     31   USE iom 
    3032   USE xios 
    3133   USE cpl_oasis3 
     
    3638   USE asmbgc, ONLY: mld_choice_bgc 
    3739   USE lbclnk 
     40   USE fabm_types, ONLY: standard_variables 
    3841 
    3942   !USE fldread         !  time interpolation 
    40  
    41    USE fabm 
    4243 
    4344   IMPLICIT NONE 
     
    4849   PRIVATE 
    4950 
    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. 
    8076 
    8177   !!---------------------------------------------------------------------- 
     
    9692      ! 
    9793      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 
    10097 
    10198!!---------------------------------------------------------------------- 
     
    109106      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 
    110107 
     108      IF (.NOT. started) CALL nemo_fabm_start 
     109 
    111110      CALL update_inputs( kt ) 
    112111 
    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 ) 
    116115 
    117116      CALL st2d_fabm_nxt( kt ) 
     
    123122      CALL trc_bc_read  ( kt )       ! tracers: surface and lateral Boundary Conditions 
    124123      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 
    125148 
    126149      IF( l_trdtrc ) THEN      ! Save the trends in the mixed layer 
     
    139162      INTEGER, INTENT(IN) :: kt 
    140163      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 
    142166       
    143167      IF (kt == nittrc000) THEN 
     
    148172      PHYT_AVG(:,:)  = 0.0 
    149173         
    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) 
    152176       
    153177      SELECT CASE( mld_choice_bgc ) 
     
    220244   END SUBROUTINE asmdiags_fabm 
    221245 
    222    SUBROUTINE compute_fabm() 
     246   SUBROUTINE compute_fabm( kt ) 
     247      INTEGER, INTENT(in) :: kt   ! ocean time-step index 
     248 
    223249      INTEGER :: ji,jj,jk,jn 
    224       TYPE(type_state) :: valid_state 
     250      LOGICAL :: valid, repaired 
    225251      REAL(wp) :: zalfg,zztmpx,zztmpy 
    226252 
    227253      ! Validate current model state (setting argument to .TRUE. enables repair=clipping) 
    228       valid_state = check_state(.TRUE.) 
    229       IF (.NOT.valid_state%valid) THEN 
     254      CALL check_state(.TRUE., valid, repaired) 
     255      IF (.NOT. valid) THEN 
    230256         WRITE(numout,*) "Invalid value in FABM encountered in area ",narea,"!!!" 
    231257#if defined key_iomput 
     
    240266#endif 
    241267      END IF 
    242       IF (valid_state%repaired) THEN 
     268      IF (repaired) THEN 
    243269         WRITE(numout,*) "Total interior repairs up to now on process",narea,":",repair_interior_count 
    244270         WRITE(numout,*) "Total surface repairs up to now on process",narea,":",repair_surface_count 
     
    246272      ENDIF 
    247273 
     274      daynumber_in_year = fjulday - fjulstartyear + 1 
     275 
    248276      ! Compute the now hydrostatic pressure 
    249277      ! copied from istate.F90 
    250278      ! ------------------------------------ 
    251279 
    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(:,:,:) 
    293314 
    294315      ! Zero rate array of interface-attached state variables 
     
    296317 
    297318      ! Compute interfacial source terms and fluxes 
    298       DO jj=1,jpj 
    299          ! Process bottom (fabm_do_bottom increments 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) 
    300321         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:)) 
    302323         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) 
    311331         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. 
    313334         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() 
    325349   END SUBROUTINE compute_fabm 
    326350 
    327    FUNCTION check_state(repair) RESULT(exit_state) 
    328       LOGICAL, INTENT(IN) :: repair 
    329       TYPE(type_state) :: exit_state 
    330  
    331       INTEGER             :: jj,jk 
    332       LOGICAL             :: valid_int,valid_sf,valid_bt 
    333  
    334       exit_state%valid = .TRUE. 
    335       exit_state%repaired =.FALSE. 
    336       DO jk=1,jpk 
    337          DO jj=1,jpj 
    338             CALL fabm_check_state(model,1,jpi,jj,jk,repair,valid_int) 
    339             IF (repair.AND..NOT.valid_int) THEN 
     351   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 
    340364               repair_interior_count = repair_interior_count + 1 
    341                exit_state%repaired = .TRUE. 
     365               repaired = .TRUE. 
    342366            END IF 
    343             IF (.NOT.(valid_int.OR.repair)) exit_state%valid = .FALSE. 
    344          END DO 
    345       END DO 
    346       DO jj=1,jpj 
    347          CALL fabm_check_surface_state(model,1,jpi,jj,repair,valid_sf) 
    348          IF (repair.AND..NOT.valid_sf) THEN 
     367            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 
    349373            repair_surface_count = repair_surface_count + 1 
    350             exit_state%repaired = .TRUE. 
     374            repaired = .TRUE. 
    351375         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) THEN 
     376         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 
    355379            repair_bottom_count = repair_bottom_count + 1 
    356             exit_state%repaired = .TRUE. 
     380            repaired = .TRUE. 
    357381         END IF 
    358          IF (.NOT.(valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE. 
    359       END DO 
    360    END FUNCTION 
     382         IF (.NOT. (valid_sf .AND. valid_bt) .AND. .NOT. repair) valid = .FALSE. 
     383      END DO 
     384   END SUBROUTINE 
    361385 
    362386   SUBROUTINE trc_sms_fabm_check_mass() 
    363387      REAL(wp) :: total(SIZE(model%conserved_quantities)) 
    364       INTEGER :: jk,jj,jn 
     388      INTEGER :: ji,jk,jj,jn 
    365389 
    366390      total = 0._wp 
    367391 
    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) 
    371397            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) 
    379407         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 
    381411         END DO 
    382412      END DO 
     
    434464 
    435465   INTEGER FUNCTION trc_sms_fabm_alloc() 
    436       INTEGER :: jj,jk,jn 
     466      INTEGER :: jn 
    437467      !!---------------------------------------------------------------------- 
    438468      !!              ***  ROUTINE trc_sms_fabm_alloc  *** 
     
    441471      ! ALLOCATE here the arrays specific to FABM 
    442472      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)) 
    446476      ! ALLOCATE( tab(...) , STAT=trc_sms_fabm_alloc ) 
    447477 
     
    452482 
    453483      ! 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)) 
    459485 
    460486      ! 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))) 
    466489#if defined key_trdtrc && defined key_iomput 
    467490      IF( lk_trdtrc ) ALLOCATE(tr_vmv(jpi,jpj,jpk,jp_fabm)) 
     
    474497      ! 
    475498 
    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 
    496515      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)) 
    498518      END DO 
    499519      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)) 
    501522      END DO 
    502523 
    503524      ! 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(:,:)) 
    523539 
    524540      ! 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.) 
    541543 
    542544      ! Set mask for negativity corrections to the relevant states 
    543       lk_rad_fabm = .FALSE. 
     545      lk_rad_fabm(:) = .FALSE. 
    544546      DO jn=1,jp_fabm 
    545         IF (model%state_variables(jn)%minimum.ge.0) THEN 
    546           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.' 
    548550        END IF 
    549551      END DO 
    550552 
    551       ! Mask land points in states with zeros, not nice, but coherent 
    552       ! with NEMO "convention": 
    553       DO jn=jp_fabm0,jp_fabm1 
    554         WHERE (tmask==0._wp) 
    555           trn(:,:,:,jn)=0._wp 
    556         END WHERE 
    557       END DO 
    558       DO jn=1,jp_fabm_surface+jp_fabm_bottom 
    559         WHERE (tmask(:,:,1)==0._wp) 
    560           fabm_st2Dn(:,:,jn)=0._wp 
    561         END WHERE 
    562       END DO 
    563553 
    564554      ! Copy initial condition for interface-attached state variables to "previous" state field 
     
    566556      fabm_st2Db = fabm_st2Dn 
    567557 
    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 
    572561 
    573562   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 
    574583 
    575584#else 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90

    r10270 r13241  
    44   !!    fabm :   Output of FABM tracers 
    55   !!====================================================================== 
    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 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_top && key_fabm && defined key_iomput 
     
    1819   USE par_fabm 
    1920   USE st2d_fabm 
    20    USE fabm, only: fabm_get_bulk_diagnostic_data, fabm_get_horizontal_diagnostic_data 
    2121 
    2222   IMPLICIT NONE 
     
    3131       MODULE PROCEDURE wri_fabm,wri_fabm_fl 
    3232   END INTERFACE trc_wri_fabm 
    33  
    3433 
    3534   PUBLIC trc_wri_fabm  
     
    5857! depth integrated 
    5958! for strict budgetting write this out at end of timestep as an average between 'now' and 'after' at kt 
    60       DO jn = 1, jp_fabm1 
    61         IF(ln_trdtrc (jn))THEN 
    62          trpool(:,:,:) = 0.5 * ( trn(:,:,:,jp_fabm0+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(:,:,:) + & 
    6362                             tr_temp(:,:,:,jn)*fse3t(:,:,:) ) 
    64          cltra = TRIM( model%state_variables(jn)%name )//"_e3t"     ! depth integrated output 
     63         cltra = TRIM( model%interior_state_variables(jn)%name )//"_e3t"     ! depth integrated output 
    6564         IF( kt == nittrc000 ) write(6,*)'output pool ',cltra 
    6665         CALL iom_put( cltra, trpool) 
     
    8079      !!--------------------------------------------------------------------- 
    8180      INTEGER, INTENT( in )               :: kt 
    82       INTEGER              :: jn 
     81      INTEGER              :: jn, jk 
     82      REAL(wp), DIMENSION(jpi,jpj)    :: vint 
    8383 
    8484#if defined key_tracer_budget 
     
    9090#endif 
    9191      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 
    93103      END DO 
    94104      DO jn = 1, jp_fabm_surface 
     
    99109      END DO 
    100110 
    101       ! write 3D diagnostics in the file 
    102       ! --------------------------------------- 
    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 DO 
    107  
    108       ! write 2D diagnostics in the file 
    109       ! --------------------------------------- 
    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 DO 
    114       ! 
    115111      CALL trc_sms_fabm_check_mass 
    116112 
     
    121117   !!  Dummy module :                                     No passive tracer 
    122118   !!---------------------------------------------------------------------- 
     119   INTERFACE trc_wri_fabm 
     120       MODULE PROCEDURE wri_fabm,wri_fabm_fl 
     121   END INTERFACE trc_wri_fabm 
     122 
    123123   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 
    127136#endif 
    128137 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90

    r12352 r13241  
    1515   USE trc 
    1616   USE par_fabm 
    17    USE fabm 
    1817   USE dom_oce 
    1918#if defined key_trdtrc && defined key_iomput 
     
    2524 
    2625#  include "domzgr_substitute.h90" 
     26#  include "vectopt_loop_substitute.h90" 
    2727 
    2828   PRIVATE 
     
    3232   ! Work arrays for vertical advection (residual movement/sinking/floating) 
    3333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: w_ct 
    34    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: w_if 
    35    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: zwgt_if 
    36    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: flux_if 
    3734#if defined key_trdtrc && defined key_iomput 
    3835   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:,:) :: tr_vmv 
     
    4138   CONTAINS 
    4239 
    43    SUBROUTINE compute_vertical_movement( kt ) 
     40   SUBROUTINE compute_vertical_movement( kt, method ) 
    4441      !!---------------------------------------------------------------------- 
    4542      !!                     ***  compute_vertical_movement  *** 
    4643      !! 
    47       !! ** Purpose :   compute vertical movement of FABM tracers 
     44      !! ** Purpose : compute vertical movement of FABM tracers through the water 
     45      !!              (sinking/floating/active movement) 
    4846      !! 
    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. 
    5449      !!---------------------------------------------------------------------- 
    5550      ! 
    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) 
    6256#if defined key_trdtrc 
    6357      CHARACTER (len=20) :: cltra 
     
    6963 
    7064      IF( neuler == 0 .AND. kt == nittrc000 ) THEN 
    71           z2dt = rdt                  ! set time step size (Euler) 
     65         z2dt = rdt                  ! set time step size (Euler) 
    7266      ELSE 
    73           z2dt = 2._wp * rdt          ! set time step size (Leapfrog) 
     67         z2dt = 2._wp * rdt          ! set time step size (Leapfrog) 
    7468      ENDIF 
     69 
    7570      ! Compute interior vertical velocities and include them in source array. 
    76       DO jj=1,jpj ! j-loop 
    77          ! 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,:)) 
    8075         END DO 
    81  
    82          DO ji=1,jpi ! i-loop 
     76         DO ji=fs_2,fs_jpim1 ! i-loop 
    8377            ! Only process this horizontal point (ji,jj) if number of layers exceeds 1 
    84             IF (mbkt(ji,jj)>1) THEN ! Level check 
    85                k_floor=mbkt(ji,jj) 
     78            k_floor = mbkt(ji,jj) 
     79            IF (k_floor > 1) THEN 
    8680               ! Linearly interpolate to velocities at the interfaces between layers 
    8781               ! Note: 
    88                !    - interface k sits between cell centre k and k-1, 
    89                !    - k [1,jpk] increases downwards 
     82               !    - interface k sits between cell centre k and k+1 (k=0 for surface) 
     83               !    - k [1,jpkm1] increases downwards 
    9084               !    - 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)) 
    10087 
    10188               ! Advect: 
    10289               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)) 
    205100                  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 
    208110            END IF ! Level check 
    209111         END DO ! i-loop 
    210112      END DO ! j-loop 
     113 
    211114#if defined key_trdtrc && defined key_iomput 
    212115      DO jn=1,jp_fabm ! State loop 
     
    220123   END SUBROUTINE compute_vertical_movement 
    221124 
     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 
    222279#endif 
    223280END MODULE 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1/NEMOGCM/NEMO/TOP_SRC/trcini.F90

    r12432 r13241  
    7676      ! +++>>> FABM 
    7777      ! Allow FABM to update numbers of biogeochemical tracers, diagnostics (jptra etc.) 
    78       IF( lk_fabm ) CALL nemo_fabm_init 
     78      IF( lk_fabm ) CALL nemo_fabm_configure 
    7979      ! FABM <<<+++ 
    8080 
     
    123123      ! Initialisation of tracers Initial Conditions 
    124124      IF( ln_trcdta )      CALL trc_dta_init(jptra) 
    125  
    126125 
    127126      IF( ln_rsttr ) THEN 
     
    162161      ! FABM +++>>> 
    163162! 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)  
    170164      ! FABM <<<+++ 
    171   
     165 
    172166      tra(:,:,:,:) = 0._wp 
    173167      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  
    325325      sn_tracer(:)%llcbc = .FALSE. 
    326326      sn_tracer(:)%llcbc = .FALSE. 
     327      sn_tracer(:)%clsname = 'NONAME' 
    327328#endif 
    328329 
     
    335336902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc in configuration namelist', lwp ) 
    336337      IF(lwm) WRITE ( numont, namtrc ) 
     338 
     339      ! +++>>> FABM 
     340      if (lk_fabm) CALL trc_nam_fabm_override(sn_tracer) 
     341      ! FABM <<<+++ 
    337342 
    338343      DO jn = 1, jptra 
     
    354359      END DO 
    355360      
    356       ! +++>>> FABM 
    357       if (lk_fabm) CALL trc_nam_fabm_override 
    358       ! FABM <<<+++ 
    359361    END SUBROUTINE trc_nam_trc 
    360362 
     
    380382      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    381383      !!--------------------------------------------------------------------- 
    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,*) '~~~~~~~' 
    386384 
    387385      IF(lwp) WRITE(numout,*) 
Note: See TracChangeset for help on using the changeset viewer.