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 13576 – NEMO

Changeset 13576


Ignore:
Timestamp:
2020-10-09T12:35:11+02:00 (3 years ago)
Author:
dford
Message:

Update NEMO-FABM coupler for FABM v1, and introduce two-way NEMO-ERSEM coupling options. See https://code.metoffice.gov.uk/trac/utils/ticket/366.

Location:
branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO
Files:
17 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90

    r13318 r13576  
    11221122#elif defined key_fabm 
    11231123      totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) + & 
    1124          &                  fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) / pnumtimes_tavg 
     1124         &                  model%get_interior_diagnostic_data(jp_fabm_o3ta) / pnumtimes_tavg 
    11251125      totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:) 
    11261126#endif 
     
    11671167         cchl_p_tavg(:,:,:)     = cchl_p(:,:,:) 
    11681168#elif defined key_fabm 
    1169          totalk_tavg(:,:,:)  = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) 
     1169         totalk_tavg(:,:,:)  = model%get_interior_diagnostic_data(jp_fabm_o3ta) 
    11701170         totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:) 
    11711171#endif 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90

    r10390 r13576  
    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         IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'25h')) THEN 
     212            fabm_3d_25h(:,:,:,jn) = model%get_interior_diagnostic_data(jn) 
     213         ENDIF 
    214214      END DO 
    215215      DO jn = 1, jp_fabm_surface 
     
    220220      END DO 
    221221      DO jn = 1, jp_fabm_2d 
    222          fabm_2d_25h(:,:,jn) = fabm_get_horizontal_diagnostic_data(model,jn) 
     222         IF (iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h')) THEN 
     223            fabm_2d_25h(:,:,jn) = model%get_horizontal_diagnostic_data(jn) 
     224         ENDIF 
    223225      END DO 
    224226#endif 
     
    327329      END DO 
    328330      DO jn = 1, jp_fabm_3d 
    329          fabm_3d_25h(:,:,:,jn) = fabm_3d_25h(:,:,:,jn) + fabm_get_interior_diagnostic_data(model, jn) 
     331         IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'25h')) THEN 
     332            fabm_3d_25h(:,:,:,jn) = fabm_3d_25h(:,:,:,jn) + model%get_interior_diagnostic_data(jn) 
     333         ENDIF 
    330334      END DO 
    331335      DO jn = 1, jp_fabm_surface 
     
    336340      END DO 
    337341      DO jn = 1, jp_fabm_2d 
    338          fabm_2d_25h(:,:,jn) = fabm_2d_25h(:,:,jn) + fabm_get_horizontal_diagnostic_data(model,jn) 
     342         IF (iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h')) THEN 
     343            fabm_2d_25h(:,:,jn) = fabm_2d_25h(:,:,jn) + model%get_horizontal_diagnostic_data(jn) 
     344         ENDIF 
    339345      END DO 
    340346#endif 
     
    401407            DO jn = 1, jp_fabm 
    402408               zw3d(:,:,:) = fabm_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    403                CALL iom_put( TRIM(model%state_variables(jn)%name)//"25h", zw3d  ) 
     409               CALL iom_put( TRIM(model%interior_state_variables(jn)%name)//"25h", zw3d  ) 
    404410            END DO 
    405411            DO jn = 1, jp_fabm_3d 
    406                zw3d(:,:,:) = fabm_3d_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    407                CALL iom_put( TRIM(model%diagnostic_variables(jn)%name)//"25h", zw3d  ) 
     412               IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'25h')) THEN 
     413                  zw3d(:,:,:) = fabm_3d_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     414                  CALL iom_put( TRIM(model%interior_diagnostic_variables(jn)%name)//"25h", zw3d  ) 
     415               ENDIF 
    408416            END DO 
    409417            DO jn = 1, jp_fabm_surface 
     
    416424            END DO 
    417425            DO jn = 1, jp_fabm_2d 
    418                zw2d(:,:) = fabm_2d_25h(:,:,jn)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
    419                CALL iom_put( TRIM(model%horizontal_diagnostic_variables(jn)%name)//"25h", zw2d  ) 
     426               IF (iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h')) THEN 
     427                  zw2d(:,:) = fabm_2d_25h(:,:,jn)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
     428                  CALL iom_put( TRIM(model%horizontal_diagnostic_variables(jn)%name)//"25h", zw2d  ) 
     429               ENDIF 
    420430            END DO 
    421431#endif 
     
    468478      END DO 
    469479      DO jn = 1, jp_fabm_3d 
    470          fabm_3d_25h(:,:,:,jn) = fabm_get_interior_diagnostic_data(model, jn) 
     480         IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'25h')) THEN 
     481            fabm_3d_25h(:,:,:,jn) = model%get_interior_diagnostic_data(jn) 
     482         ENDIF 
    471483      END DO 
    472484      DO jn = 1, jp_fabm_surface 
     
    477489      END DO 
    478490      DO jn = 1, jp_fabm_2d 
    479          fabm_2d_25h(:,:,jn) = fabm_get_horizontal_diagnostic_data(model,jn) 
     491         IF (iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h')) THEN 
     492            fabm_2d_25h(:,:,jn) = model%get_horizontal_diagnostic_data(jn) 
     493         ENDIF 
    480494      END DO 
    481495#endif 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90

    r10390 r13576  
    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            IF ( iom_use('top_'//TRIM(model%interior_diagnostic_variables(jn)%name)) .OR. & 
     179               & iom_use('mid_'//TRIM(model%interior_diagnostic_variables(jn)%name)) .OR. & 
     180               & iom_use('bot_'//TRIM(model%interior_diagnostic_variables(jn)%name)) ) THEN 
     181               CALL dia_calctmb( model%get_interior_diagnostic_data(jn), zwtmb ) 
     182               CALL iom_put( "top_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,1) ) 
     183               CALL iom_put( "mid_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,2) ) 
     184               CALL iom_put( "bot_"//TRIM(model%interior_diagnostic_variables(jn)%name) , zwtmb(:,:,3) ) 
     185            ENDIF 
    183186         END DO 
    184187#endif 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r8059 r13576  
    3232   USE wrk_nemo       ! Memory Allocation 
    3333   USE timing         ! Timing 
     34#if defined key_fabm 
     35   USE trc, ONLY: trn  ! FABM variables 
     36   USE par_fabm        ! FABM parameters 
     37#endif 
    3438 
    3539   IMPLICIT NONE 
     
    4549   LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag 
    4650   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem) 
     51   LOGICAL , PUBLIC ::   ln_qsr_spec  !: spectral model heating from ERSEM 
    4752   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0) 
    4853   INTEGER , PUBLIC ::   nn_kd490dta  !: use kd490dta data (=1) or not (=0) 
     
    106111      REAL(wp) ::   zz0, zz1, z1_e3t     !    -         - 
    107112      REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr 
     113      REAL(wp), POINTER, DIMENSION(:,:,:) :: zekb_3d, zekg_3d, zekr_3d 
    108114      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt 
    109115      !!---------------------------------------------------------------------- 
     
    113119      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )  
    114120      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
     121      CALL wrk_alloc( jpi, jpj, jpk, zekb_3d, zekg_3d, zekr_3d )   
    115122      ! 
    116123      IF( kt == nit000 ) THEN 
     
    150157       
    151158      !                                           ! ============================================== ! 
    152       IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  ! 
     159      IF( ln_qsr_spec ) THEN                      !  ERSEM spectral heating                        ! 
     160         !                                        ! ============================================== ! 
     161         DO jk = 1, jpkm1 
     162           qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) ) 
     163         END DO 
     164         !                                        Add to the general trend 
     165         DO jk = 1, jpkm1 
     166            DO jj = 2, jpjm1 
     167               DO ji = fs_2, fs_jpim1   ! vector opt. 
     168                  z1_e3t = zfact / fse3t(ji,jj,jk) 
     169                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t 
     170               END DO 
     171            END DO 
     172         END DO 
     173         zea(:,:,:) = etot3(:,:,:) * tmask(:,:,:) 
     174         DO jk = jpkm1, 1, -1 
     175            zea(:,:,jk) = zea(:,:,jk) + zea(:,:,jk+1) 
     176         END DO 
     177         CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution 
     178         IF ( ln_qsr_ice ) THEN 
     179            DO jj = 1, jpj 
     180               DO ji = 1, jpi 
     181                  IF ( qsr(ji,jj) /= 0._wp ) THEN 
     182                     fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) ) 
     183                  ELSE 
     184                     fraqsr_1lev(ji,jj) = 1. 
     185                  ENDIF 
     186               END DO 
     187            END DO 
     188         ENDIF 
     189         !  
     190 
     191         !                                        ! ============================================== ! 
     192      ELSEIF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  ! 
    153193         !                                        ! ============================================== ! 
    154194         DO jk = 1, jpkm1 
     
    185225            !                                             ! ------------------------- ! 
    186226            ! Set chlorophyl concentration 
    187             IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume 
    188                ! 
    189                IF( nn_chldta == 1 ) THEN                             !*  Variable Chlorophyll 
     227            IF( nn_chldta == 2 .OR. nn_chldta == 1 .OR. lk_vvl ) THEN   !*  Variable Chlorophyll or ocean volume 
     228               ! 
     229               IF( nn_chldta == 2 ) THEN 
     230                  DO jk = 1, nksr+1 
     231                     DO jj = 1, jpj 
     232                        DO ji = 1, jpi 
     233#if defined key_fabm 
     234                           zchl = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) + & 
     235                              &   trn(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) 
     236#endif 
     237                           zchl = MIN(  10. , MAX( 0.03, zchl )  ) 
     238                           irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
     239                           !                                                          
     240                           zekb_3d(ji,jj,jk) = rkrgb(1,irgb) 
     241                           zekg_3d(ji,jj,jk) = rkrgb(2,irgb) 
     242                           zekr_3d(ji,jj,jk) = rkrgb(3,irgb) 
     243                        END DO 
     244                     END DO 
     245                  END DO 
     246                  ! 
     247               ELSEIF( nn_chldta == 1 ) THEN                             !*  Variable Chlorophyll 
    190248                  ! 
    191249                  CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step 
     
    219277               ! 
    220278               DO jk = 2, nksr+1 
     279                  IF( nn_chldta == 2 ) THEN 
     280                     zekb(:,:) = zekb_3d(:,:,jk) 
     281                     zekg(:,:) = zekg_3d(:,:,jk) 
     282                     zekr(:,:) = zekr_3d(:,:,jk) 
     283                  ENDIF 
    221284!CDIR NOVERRCHK 
    222285                  DO jj = 1, jpj 
     
    237300               ! clem: store attenuation coefficient of the first ocean level 
    238301               IF ( ln_qsr_ice ) THEN 
     302                  IF( nn_chldta == 2 ) THEN 
     303                     zekb(:,:) = zekb_3d(:,:,1) 
     304                     zekg(:,:) = zekg_3d(:,:,1) 
     305                     zekr(:,:) = zekr_3d(:,:,1) 
     306                  ENDIF 
    239307                  DO jj = 1, jpj 
    240308                     DO ji = 1, jpi 
     
    386454      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )  
    387455      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
     456      CALL wrk_dealloc( jpi, jpj, jpk, zekb_3d, zekg_3d, zekr_3d )  
    388457      ! 
    389458      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr') 
     
    423492      !! 
    424493      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 
     494         &                  ln_qsr_spec, nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta 
    426495      !!---------------------------------------------------------------------- 
    427496 
     
    451520         WRITE(numout,*) '      2 band               light penetration   ln_qsr_2bd = ', ln_qsr_2bd 
    452521         WRITE(numout,*) '      bio-model            light penetration   ln_qsr_bio = ', ln_qsr_bio 
     522         WRITE(numout,*) '      ERSEM spectral heating model             ln_qsr_spec= ', ln_qsr_spec 
    453523         WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice 
    454          WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)    nn_chldta = ', nn_chldta 
     524         WRITE(numout,*) '      RGB: model (2), file (1) or cst (0) chl   nn_chldta = ', nn_chldta 
    455525         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs = ', rn_abs 
    456526         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0 
     
    470540         IF( ln_qsr_2bd  )   ioptio = ioptio + 1 
    471541         IF( ln_qsr_bio  )   ioptio = ioptio + 1 
     542         IF( ln_qsr_spec )   ioptio = ioptio + 1 
    472543         IF( nn_kd490dta == 1 )   ioptio = ioptio + 1 
    473544         ! 
     
    478549         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1  
    479550         IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2 
    480          IF( ln_qsr_2bd                      )   nqsr =  3 
    481          IF( ln_qsr_bio                      )   nqsr =  4 
    482          IF( nn_kd490dta == 1                )   nqsr =  5 
     551         IF( ln_qsr_rgb .AND. nn_chldta == 2 )   nqsr =  3 
     552         IF( ln_qsr_2bd                      )   nqsr =  4 
     553         IF( ln_qsr_bio                      )   nqsr =  5 
     554         IF( nn_kd490dta == 1                )   nqsr =  6 
     555         IF( ln_qsr_spec                     )   nqsr =  7 
    483556         ! 
    484557         IF(lwp) THEN                   ! Print the choice 
    485558            WRITE(numout,*) 
    486559            IF( nqsr ==  1 )   WRITE(numout,*) '         R-G-B   light penetration - Constant Chlorophyll' 
    487             IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B   light penetration - Chl data ' 
    488             IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration' 
    489             IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration' 
    490             IF( nqsr ==  5 )   WRITE(numout,*) '         KD490 light penetration' 
    491          ENDIF 
     560            IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B   light penetration - Chl data from file' 
     561            IF( nqsr ==  3 )   WRITE(numout,*) '         R-G-B   light penetration - Chl data from model' 
     562            IF( nqsr ==  4 )   WRITE(numout,*) '         2 bands light penetration' 
     563            IF( nqsr ==  5 )   WRITE(numout,*) '         bio-model light penetration' 
     564            IF( nqsr ==  6 )   WRITE(numout,*) '         KD490 light penetration' 
     565            IF( nqsr ==  7 )   WRITE(numout,*) '         ERSEM spectral light penetration' 
     566         ENDIF 
     567#if ! defined key_fabm 
     568         ! 
     569         IF( nqsr ==  2 ) THEN 
     570            CALL ctl_stop( 'nn_chldta=2 so trying to use ERSEM chlorophyll for light penetration', & 
     571               &           'but not running with ERSEM' ) 
     572         ELSEIF( nqsr ==  7 ) THEN 
     573            CALL ctl_stop( 'ln_qsr_spec=.true. so trying to use ERSEM spectral light penetration', & 
     574               &           'but not running with ERSEM' ) 
     575         ENDIF 
     576#endif 
    492577         ! 
    493578      ENDIF 
     
    536621                  &                                         'Solar penetration function of read chlorophyll', 'namtra_qsr' ) 
    537622               ! 
     623            ELSEIF( nn_chldta == 2 ) THEN       !* Chl data will be got from model at each time step 
     624               IF(lwp) WRITE(numout,*) 
     625               IF(lwp) WRITE(numout,*) '        Chlorophyll will be taken from model at each time step' 
    538626            ELSE                                !* constant Chl : compute once for all the distribution of light (etot3) 
    539627               IF(lwp) WRITE(numout,*) 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r11277 r13576  
    158158         END DO 
    159159#else 
    160           IF( lk_asminc ) THEN 
    161              IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 )    ! Output background fields 
    162              IF( ln_asmdin ) THEN                        ! Direct initialization 
    163                 IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 )    ! Tracers 
    164                 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics 
    165                 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH 
    166                 IF( lk_bgcinc ) CALL bgc_asm_inc( nit000 - 1 )    ! BGC 
    167              ENDIF 
    168           ENDIF 
     160          ! Initial call to assimilation routines moved to stp 
    169161 
    170162#if defined key_agrif 
     
    493485      IF( lk_diaobs     ) THEN                  ! Observation & model comparison 
    494486                            CALL dia_obs_init            ! Initialize observational data 
    495                             CALL dia_obs( nit000 - 1 )   ! Observation operator for restart 
    496       ENDIF 
    497  
    498       !                                     ! Assimilation increments 
    499       IF( lk_asminc ) THEN  
    500 #if defined key_shelf  
    501          CALL  zdf_mxl_tref()     ! Initialization of hmld_tref 
    502 #endif  
    503          CALL asm_inc_init     ! Initialize assimilation increments  
    504       ENDIF 
     487                            ! Initial call to dia_obs moved to stp 
     488      ENDIF 
     489 
     490      ! Initialisation of assimilation moved to stp 
    505491 
    506492      IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler 
    507                             CALL dia_tmb_init  ! TMB outputs 
    508                             CALL dia_25h_init  ! 25h mean  outputs 
    509                             CALL dia_diaopfoam_init  ! FOAM operational output 
     493       
     494      ! Initialisation of tmb/25h/diaopfoam outputs moved to stp 
     495       
    510496      ! 
    511497   END SUBROUTINE nemo_init 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/step.F90

    r11277 r13576  
    9696                      CALL iom_init(      cxios_context          )  ! iom_put initialization 
    9797         IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" )  ! initialize context for coarse grid 
     98          
     99#if defined key_fabm 
     100         ! FABM can only finish initialising once IOM has 
     101         IF ( lk_fabm ) CALL nemo_fabm_start 
     102#endif 
     103          
     104         ! First call to dia_* and asm_inc_init must wait for FABM to be initialised 
     105         ! (if running with FABM, but no harm moving the calls to here from nemo_init either way) 
     106         CALL dia_obs( nit000 - 1 )   ! Observation operator for restart 
     107          
     108         IF( lk_asminc ) THEN  
     109#if defined key_shelf 
     110            CALL  zdf_mxl_tref()     ! Initialization of hmld_tref 
     111#endif  
     112            CALL asm_inc_init     ! Initialize assimilation increments  
     113            IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 )    ! Output background fields 
     114            IF( ln_asmdin ) THEN                        ! Direct initialization 
     115               IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 )    ! Tracers 
     116               IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics 
     117               IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH 
     118               IF( lk_bgcinc ) CALL bgc_asm_inc( nit000 - 1 )    ! BGC 
     119            ENDIF 
     120         ENDIF 
     121          
     122         CALL dia_tmb_init        ! TMB outputs 
     123         CALL dia_25h_init        ! 25h mean  outputs 
     124         CALL dia_diaopfoam_init  ! FOAM operational output 
    98125      ENDIF 
    99126 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r11277 r13576  
    126126   USE trcstp           ! passive tracer time-stepping      (trc_stp routine) 
    127127#endif 
     128#if defined key_fabm 
     129   USE par_fabm, ONLY: &    ! FABM parameters 
     130      & lk_fabm 
     131   USE trcsms_fabm, ONLY: & ! FABM routines 
     132      & nemo_fabm_start 
     133#endif 
     134   USE diatmb          ! Top,middle,bottom output 
     135   USE dia25h          ! 25h mean output 
     136   USE diaopfoam       ! FOAM operational output 
    128137   !!---------------------------------------------------------------------- 
    129138   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/inputs_fabm.F90

    r10158 r13576  
    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/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90

    r11545 r13576  
    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, & 
     
    3539                      jp_fabm_r8c,   jp_fabm_r8p,  & 
    3640                      jp_fabm_r8s,   jp_fabm_xeps, & 
     41                      jp_fabm_swr,   jp_fabm_kd490, & 
    3742                      jp_fabm_pgrow, jp_fabm_ploss 
    3843 
     44   LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) ::   lk_rad_fabm !: FABM negativity correction flag array 
     45 
    3946#if defined key_fabm 
     47   CLASS (type_fabm_model), POINTER :: model !FABM model instance 
     48 
    4049   !!--------------------------------------------------------------------- 
    4150   !!   'key_fabm'                     FABM tracers 
    4251   !!--------------------------------------------------------------------- 
    4352   LOGICAL, PUBLIC, PARAMETER ::   lk_fabm     = .TRUE.   !: FABM flag  
    44    LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) ::   lk_rad_fabm !: FABM negativity correction flag array  
    4553#else 
    4654   !!--------------------------------------------------------------------- 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90

    r12352 r13576  
    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 
     
    123134      jp_fabm_o3pc  = fabm_diag_index( 'O3_pCO2' ) 
    124135      jp_fabm_xeps  = fabm_diag_index( 'light_xEPS' ) 
     136      jp_fabm_swr   = fabm_diag_index( 'light_swr_abs' ) 
     137      jp_fabm_kd490 = fabm_diag_index( 'light_Kd_band3' ) 
    125138      jp_fabm_pgrow = fabm_diag_index( 'p_grow_sum_result' ) 
    126139      jp_fabm_ploss = fabm_diag_index( 'p_loss_sum_result' ) 
     
    128141      IF (lwp) THEN 
    129142         ! write field_def_fabm.xml on lead process 
    130          OPEN(UNIT=xml_unit,FILE='field_def_fabm.xml',ACTION='WRITE',STATUS='REPLACE') 
     143         OPEN(UNIT=xml_unit, FILE='field_def_fabm.xml', ACTION='WRITE', STATUS='REPLACE') 
    131144 
    132145         WRITE (xml_unit,1000) '<field_definition level="1" prec="4" operation="average" enabled=".TRUE." default_value="1.e20" >' 
     
    134147         WRITE (xml_unit,1000) ' <field_group id="ptrc_T" grid_ref="grid_T_3D">' 
    135148         DO jn=1,jp_fabm 
    136             CALL write_variable_xml(xml_unit,model%state_variables(jn)) 
     149            CALL write_variable_xml(xml_unit,model%interior_state_variables(jn)) 
    137150#if defined key_trdtrc 
    138             CALL write_trends_xml(xml_unit,model%state_variables(jn)) 
    139 #endif 
    140             CALL write_25hourm_xml(xml_unit,model%state_variables(jn)) 
    141             CALL write_tmb_xml(xml_unit,model%state_variables(jn)) 
     151            CALL write_trends_xml(xml_unit,model%interior_state_variables(jn)) 
     152#endif 
     153            CALL write_25hourm_xml(xml_unit,model%interior_state_variables(jn)) 
     154            CALL write_tmb_xml(xml_unit,model%interior_state_variables(jn)) 
    142155         END DO 
    143156         WRITE (xml_unit,1000) ' </field_group>' 
     
    155168 
    156169         WRITE (xml_unit,1000) ' <field_group id="diad_T" grid_ref="grid_T_2D">' 
    157          DO jn=1,size(model%diagnostic_variables) 
    158             CALL write_variable_xml(xml_unit,model%diagnostic_variables(jn),3) 
    159             CALL write_25hourm_xml(xml_unit,model%diagnostic_variables(jn),3) 
    160             CALL write_tmb_xml(xml_unit,model%diagnostic_variables(jn)) 
     170         DO jn=1,size(model%interior_diagnostic_variables) 
     171            CALL write_variable_xml(xml_unit,model%interior_diagnostic_variables(jn),3) 
     172            CALL write_25hourm_xml(xml_unit,model%interior_diagnostic_variables(jn),3) 
     173            CALL write_tmb_xml(xml_unit,model%interior_diagnostic_variables(jn)) 
    161174         END DO 
    162175         DO jn=1,size(model%horizontal_diagnostic_variables) 
    163176            CALL write_variable_xml(xml_unit,model%horizontal_diagnostic_variables(jn)) 
    164177            CALL write_25hourm_xml(xml_unit,model%horizontal_diagnostic_variables(jn)) 
     178         END DO 
     179         DO jn=1,size(model%interior_state_variables) 
     180            WRITE (xml_unit,'(A)') '  <field id="'//TRIM(model%interior_state_variables(jn)%name)// & 
     181               &                   '_VINT" long_name="depth-integrated '//TRIM(model%interior_state_variables(jn)%long_name)// & 
     182               &                   '" unit="'//TRIM(model%interior_state_variables(jn)%units)//'*m" default_value="0.0" />' 
     183         END DO 
     184         DO jn=1,size(model%interior_diagnostic_variables) 
     185            WRITE (xml_unit,'(A)') '  <field id="'//TRIM(model%interior_diagnostic_variables(jn)%name)// & 
     186               &                   '_VINT" long_name="depth-integrated '//TRIM(model%interior_diagnostic_variables(jn)%long_name)// & 
     187               &                   '" unit="'//TRIM(model%interior_diagnostic_variables(jn)%units)//'*m" default_value="0.0" />' 
    165188         END DO 
    166189         WRITE (xml_unit,1000) ' </field_group>' 
     
    1952181000 FORMAT (A) 
    196219 
    197    END SUBROUTINE nemo_fabm_init 
     220   END SUBROUTINE nemo_fabm_configure 
    198221 
    199222   SUBROUTINE write_variable_xml(xml_unit,variable,flag_grid_ref) 
    200223      INTEGER,INTENT(IN) :: xml_unit 
    201224      INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 
    202       CLASS (type_external_variable),INTENT(IN) :: variable 
     225      CLASS (type_fabm_variable),INTENT(IN) :: variable 
    203226 
    204227      CHARACTER(LEN=20) :: missing_value,string_dimensions 
     
    233256      INTEGER,INTENT(IN) :: xml_unit 
    234257      INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 
    235       CLASS (type_external_variable),INTENT(IN) :: variable 
     258      CLASS (type_fabm_variable),INTENT(IN) :: variable 
    236259 
    237260      CHARACTER(LEN=20) :: missing_value,string_dimensions 
     
    265288   SUBROUTINE write_tmb_xml(xml_unit,variable) 
    266289      INTEGER,INTENT(IN) :: xml_unit 
    267       CLASS (type_external_variable),INTENT(IN) :: variable 
     290      CLASS (type_fabm_variable),INTENT(IN) :: variable 
    268291 
    269292      CHARACTER(LEN=20) :: missing_value 
     
    279302      INTEGER,INTENT(IN) :: xml_unit 
    280303      INTEGER,INTENT(IN),OPTIONAL :: flag_grid_ref 
    281       CLASS (type_external_variable),INTENT(IN) :: variable 
     304      CLASS (type_fabm_variable),INTENT(IN) :: variable 
    282305 
    283306      INTEGER :: number_dimensions,i 
     
    383406      !! ** Purpose :   initialization for FABM model 
    384407      !! 
    385       !! ** Method  : - Read the namcfc namelist and check the parameter values 
     408      !! ** Method  : - Allocate FABM arrays, configure domain, send data 
    386409      !!---------------------------------------------------------------------- 
    387410#if defined key_git_version 
    388       TYPE (type_version),POINTER :: version 
     411      TYPE (type_version), POINTER :: version 
    389412#endif 
    390413      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' ) 
    394414 
    395415      IF(lwp) WRITE(numout,*) 
     
    399419      IF(lwp) WRITE(numout,*) ' NEMO version:   ',git_commit_id,' (',git_branch_name,' branch)' 
    400420      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 
    405421      version => first_module_version 
    406  
    407422      do while (associated(version)) 
    408423         IF(lwp) WRITE(numout,*)  ' '//trim(version%module_name)//' version:   ',trim(version%version_string) 
     
    411426#endif 
    412427 
     428      ! Allocate FABM arrays 
     429      IF(trc_sms_fabm_alloc() /= 0) CALL ctl_stop( 'STOP', 'trc_ini_fabm: unable to allocate FABM arrays' ) 
     430 
    413431      ! Log mapping of FABM states: 
    414432      IF (lwp) THEN 
    415          IF (jp_fabm.gt.0) WRITE(numout,*) " FABM tracers:" 
     433         IF (jp_fabm > 0) WRITE(numout,*) " FABM tracers:" 
    416434         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:" 
     435            WRITE(numout,*) "   State",jn,":",trim(model%interior_state_variables(jn)%name), & 
     436               " (",trim(model%interior_state_variables(jn)%long_name), & 
     437               ") [",trim(model%interior_state_variables(jn)%units),"]" 
     438         END DO 
     439         IF (jp_fabm_surface > 0) WRITE(numout,*) "FABM seasurface states:" 
    422440         DO jn=1,jp_fabm_surface 
    423441            WRITE(numout,*) "   State",jn,":",trim(model%surface_state_variables(jn)%name), & 
    424442               " (",trim(model%surface_state_variables(jn)%long_name), & 
    425443               ") [",trim(model%surface_state_variables(jn)%units),"]" 
    426          ENDDO 
    427          IF (jp_fabm_bottom.gt.0) WRITE(numout,*) "FABM seafloor states:" 
     444         END DO 
     445         IF (jp_fabm_bottom > 0) WRITE(numout,*) "FABM seafloor states:" 
    428446         DO jn=1,jp_fabm_bottom 
    429447            WRITE(numout,*) "   State",jn,":",trim(model%bottom_state_variables(jn)%name), & 
    430448               " (",trim(model%bottom_state_variables(jn)%long_name), & 
    431449               ") [",trim(model%bottom_state_variables(jn)%units),"]" 
    432          ENDDO 
    433       ENDIF 
     450         END DO 
     451      END IF 
    434452       
    435453      ! Initialise variables required for Hemmings et al. (2008) 
     
    442460   END SUBROUTINE trc_ini_fabm 
    443461 
     462   SUBROUTINE nemo_fabm_driver_fatal_error(self, location, message) 
     463      CLASS (type_nemo_fabm_driver), INTENT(INOUT) :: self 
     464      CHARACTER(len=*),              INTENT(IN)    :: location, message 
     465 
     466      CALL ctl_stop('STOP', TRIM(location)//': '//TRIM(message)) 
     467      STOP 
     468   END SUBROUTINE 
     469 
     470   SUBROUTINE nemo_fabm_driver_log_message(self, message) 
     471      CLASS (type_nemo_fabm_driver), INTENT(INOUT) :: self 
     472      CHARACTER(len=*),              INTENT(IN)    :: message 
     473 
     474      IF(lwp) WRITE (numout,*) TRIM(message) 
     475   END SUBROUTINE 
     476 
    444477   INTEGER FUNCTION fabm_state_index( state_name ) 
    445478      !!---------------------------------------------------------------------- 
     
    453486      IMPLICIT NONE 
    454487       
    455       CHARACTER(LEN=256), INTENT(IN) :: state_name 
    456        
    457       INTEGER                        :: jn 
     488      CHARACTER(LEN=*), INTENT(IN) :: state_name 
     489       
     490      INTEGER                      :: jn 
    458491 
    459492      !!---------------------------------------------------------------------- 
     
    461494      fabm_state_index = -1 
    462495      DO jn=1,jp_fabm 
    463          IF (TRIM(model%state_variables(jn)%name) == TRIM(state_name)) THEN 
     496         IF (TRIM(model%interior_state_variables(jn)%name) == TRIM(state_name)) THEN 
    464497            fabm_state_index = jn 
    465498            EXIT 
     
    485518      IMPLICIT NONE 
    486519       
    487       CHARACTER(LEN=256), INTENT(IN) :: diag_name 
    488        
    489       INTEGER                        :: jn 
     520      CHARACTER(LEN=*), INTENT(IN) :: diag_name 
     521       
     522      INTEGER                      :: jn 
    490523 
    491524      !!---------------------------------------------------------------------- 
    492525       
    493526      fabm_diag_index = -1 
    494       DO jn = 1, SIZE(model%diagnostic_variables) 
    495          IF (TRIM(model%diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN 
     527      DO jn = 1, SIZE(model%interior_diagnostic_variables) 
     528         IF (TRIM(model%interior_diagnostic_variables(jn)%name) == TRIM(diag_name)) THEN 
    496529            fabm_diag_index = jn 
    497530            EXIT 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcnam_fabm.F90

    r10270 r13576  
    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/NEMOGCM/NEMO/TOP_SRC/FABM/trcrst_fabm.F90

    r10156 r13576  
    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/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90

    r10728 r13576  
    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 
     41   USE par_fabm, ONLY: jp_fabm_swr 
    3842 
    3943   !USE fldread         !  time interpolation 
    40  
    41    USE fabm 
    4244 
    4345   IMPLICIT NONE 
     
    4850   PRIVATE 
    4951 
    50    PUBLIC   trc_sms_fabm       ! called by trcsms.F90 module 
    51    PUBLIC   trc_sms_fabm_alloc ! called by trcini_fabm.F90 module 
    52    PUBLIC   trc_sms_fabm_check_mass 
    53    PUBLIC   st2d_fabm_nxt ! 2D state intergration 
    54    PUBLIC   compute_fabm ! Compute FABM sources, sinks and diagnostics 
    55  
    56    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: flux    ! Cross-interface flux of pelagic variables (# m-2 s-1) 
    57  
    58    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   :: ext     ! Light extinction coefficient (m-1) 
    59  
    60    ! Work array for mass aggregation 
    61    REAL(wp), ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: current_total 
    62  
    63  
    64    ! Arrays for environmental variables 
    65    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: prn,rho 
    66    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: taubot 
    67  
    68    ! repair counters 
    69    INTEGER :: repair_interior_count,repair_surface_count,repair_bottom_count 
    70  
    71    ! state check type 
    72    TYPE type_state 
    73       LOGICAL             :: valid 
    74       LOGICAL             :: repaired 
    75    END TYPE 
    76  
    77    REAL(wp), PUBLIC :: daynumber_in_year 
    78  
    79    TYPE (type_bulk_variable_id),SAVE :: swr_id 
     52   PUBLIC   trc_sms_fabm            ! called by trcsms.F90 module 
     53   PUBLIC   trc_sms_fabm_alloc      ! called by trcini_fabm.F90 module 
     54   PUBLIC   trc_sms_fabm_check_mass ! called by trcwri_fabm.F90 
     55   PUBLIC   nemo_fabm_start 
     56 
     57   ! Work arrays 
     58   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: flux          ! Cross-interface flux of pelagic variables (# m-2 s-1) 
     59   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: current_total ! Totals of conserved quantities 
     60 
     61   ! Arrays for environmental variables that are computed by the coupler 
     62   REAL(wp), PUBLIC, TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: prn,rho 
     63   REAL(wp), PUBLIC, TARGET, ALLOCATABLE, DIMENSION(:,:) :: taubot 
     64   REAL(wp), PUBLIC, TARGET :: daynumber_in_year 
     65 
     66   ! State repair counters 
     67   INTEGER, SAVE :: repair_interior_count = 0 
     68   INTEGER, SAVE :: repair_surface_count  = 0 
     69   INTEGER, SAVE :: repair_bottom_count   = 0 
     70 
     71   ! Coupler parameters 
     72   INTEGER, PUBLIC :: nn_adv  ! Vertical advection scheme for sinking/floating/movement 
     73                              ! (1: 1st order upwind, 3: 3rd order TVD) 
     74 
     75   ! Flag indicating whether model%start has been called (will be done on-demand) 
     76   LOGICAL, SAVE :: started = .false. 
    8077 
    8178   !!---------------------------------------------------------------------- 
     
    9693      ! 
    9794      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    98       INTEGER :: jn 
    99       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrfabm 
     95      INTEGER :: jn, jk 
     96      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrfabm, pdat 
     97      REAL(wp), DIMENSION(jpi,jpj)    :: vint 
    10098 
    10199!!---------------------------------------------------------------------- 
     
    109107      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 
    110108 
     109      IF (.NOT. started) CALL nemo_fabm_start 
     110 
    111111      CALL update_inputs( kt ) 
    112112 
    113       CALL compute_fabm 
    114  
    115       CALL compute_vertical_movement( kt ) 
     113      CALL compute_fabm( kt ) 
     114 
     115      CALL compute_vertical_movement( kt, nn_adv ) 
    116116 
    117117      CALL st2d_fabm_nxt( kt ) 
     
    123123      CALL trc_bc_read  ( kt )       ! tracers: surface and lateral Boundary Conditions 
    124124      CALL trc_rnf_fabm ( kt ) ! River forcings 
     125 
     126      ! Send 3D diagnostics to output (these apply to time "n") 
     127      DO jn = 1, size(model%interior_diagnostic_variables) 
     128         IF (model%interior_diagnostic_variables(jn)%save) THEN 
     129            ! Save 3D field 
     130            pdat => model%get_interior_diagnostic_data(jn) 
     131            CALL iom_put(model%interior_diagnostic_variables(jn)%name, pdat) 
     132 
     133            ! Save depth integral if selected for output in XIOS 
     134            IF (iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT')) THEN 
     135               vint = 0._wp 
     136               DO jk = 1, jpkm1 
     137                  vint = vint + pdat(:,:,jk) * fse3t(:,:,jk) * tmask(:,:,jk) 
     138               END DO 
     139               CALL iom_put(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT', vint) 
     140            END IF 
     141         END IF 
     142      END DO 
     143 
     144      ! Send 2D diagnostics to output (these apply to time "n") 
     145      DO jn = 1, size(model%horizontal_diagnostic_variables) 
     146         IF (model%horizontal_diagnostic_variables(jn)%save) & 
     147             CALL iom_put( model%horizontal_diagnostic_variables(jn)%name, model%get_horizontal_diagnostic_data(jn)) 
     148      END DO 
    125149 
    126150      IF( l_trdtrc ) THEN      ! Save the trends in the mixed layer 
     
    139163      INTEGER, INTENT(IN) :: kt 
    140164      INTEGER :: ji,jj,jk,jkmax 
    141       REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d, zmld 
     165      REAL(wp), DIMENSION(jpi,jpj) :: zmld 
     166      REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d 
    142167       
    143168      IF (kt == nittrc000) THEN 
     
    148173      PHYT_AVG(:,:)  = 0.0 
    149174         
    150       pgrow_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_pgrow) 
    151       ploss_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_ploss) 
     175      pgrow_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_pgrow) 
     176      ploss_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_ploss) 
    152177       
    153178      SELECT CASE( mld_choice_bgc ) 
     
    220245   END SUBROUTINE asmdiags_fabm 
    221246 
    222    SUBROUTINE compute_fabm() 
     247   SUBROUTINE compute_fabm( kt ) 
     248      INTEGER, INTENT(in) :: kt   ! ocean time-step index 
     249 
    223250      INTEGER :: ji,jj,jk,jn 
    224       TYPE(type_state) :: valid_state 
     251      LOGICAL :: valid, repaired 
    225252      REAL(wp) :: zalfg,zztmpx,zztmpy 
    226253 
    227254      ! Validate current model state (setting argument to .TRUE. enables repair=clipping) 
    228       valid_state = check_state(.TRUE.) 
    229       IF (.NOT.valid_state%valid) THEN 
     255      CALL check_state(.TRUE., valid, repaired) 
     256      IF (.NOT. valid) THEN 
    230257         WRITE(numout,*) "Invalid value in FABM encountered in area ",narea,"!!!" 
    231258#if defined key_iomput 
     
    240267#endif 
    241268      END IF 
    242       IF (valid_state%repaired) THEN 
     269      IF (repaired) THEN 
    243270         WRITE(numout,*) "Total interior repairs up to now on process",narea,":",repair_interior_count 
    244271         WRITE(numout,*) "Total surface repairs up to now on process",narea,":",repair_surface_count 
     
    246273      ENDIF 
    247274 
     275      daynumber_in_year = fjulday - fjulstartyear + 1 
     276 
    248277      ! Compute the now hydrostatic pressure 
    249278      ! copied from istate.F90 
    250279      ! ------------------------------------ 
    251280 
    252       zalfg = 0.5e-4 * grav ! FABM wants dbar, convert from Pa 
    253  
    254       rho = rau0 * ( 1. + rhd ) 
    255  
    256       prn(:,:,1) = 10.1325 + zalfg * fse3t(:,:,1) * rho(:,:,1) 
    257  
    258       daynumber_in_year=(fjulday-fjulstartyear+1)*1._wp 
    259  
    260       DO jk = 2, jpk                                              ! Vertical integration from the surface 
    261          prn(:,:,jk) = prn(:,:,jk-1) + zalfg * ( & 
    262                      fse3t(:,:,jk-1) * rho(:,:,jk-1)  & 
    263                      + fse3t(:,:,jk) * rho(:,:,jk) ) 
    264       END DO 
    265  
    266       ! Bottom stress 
    267       taubot(:,:) = 0._wp 
    268       DO jj = 2, jpjm1 
    269          DO ji = fs_2, fs_jpim1   ! vector opt. 
    270                zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  & 
    271                       &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  ) 
    272                zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  & 
    273                       &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
    274                taubot(ji,jj) = 0.5_wp * rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
    275                ! 
    276          ENDDO 
    277       ENDDO 
    278       ! Compute light extinction 
    279       DO jk=1,jpk 
    280           DO jj=1,jpj 
    281             call fabm_get_light_extinction(model,1,jpi,jj,jk,ext) 
    282          END DO 
    283       END DO 
    284  
    285       ! Compute light field (stored among FABM's internal diagnostics) 
    286       DO jj=1,jpj 
    287           DO ji=1,jpi 
    288             call fabm_get_light(model,1,jpk,ji,jj) 
    289          END DO 
    290       END DO 
    291  
    292       ! TODO: retrieve 3D shortwave and store in etot3 
     281      IF (ALLOCATED(rho)) rho = rau0 * ( 1._wp + rhd ) 
     282 
     283      IF (ALLOCATED(prn)) THEN 
     284         zalfg = 0.5e-4_wp * grav ! FABM wants dbar, convert from Pa (and multiply with 0.5 to average 2 cell thicknesses below) 
     285         prn(:,:,1) = 10.1325_wp + zalfg * fse3t(:,:,1) * rho(:,:,1) 
     286         DO jk = 2, jpkm1                                              ! Vertical integration from the surface 
     287            prn(:,:,jk) = prn(:,:,jk-1) + zalfg * ( & 
     288                        fse3t(:,:,jk-1) * rho(:,:,jk-1)  & 
     289                        + fse3t(:,:,jk) * rho(:,:,jk) ) 
     290         END DO 
     291      END IF 
     292 
     293      ! Compute the bottom stress 
     294      ! copied from diawri.F90 
     295      ! ------------------------------------ 
     296 
     297      IF (ALLOCATED(taubot)) THEN 
     298         taubot(:,:) = 0._wp 
     299         DO jj = 2, jpjm1 
     300            DO ji = fs_2, fs_jpim1   ! vector opt. 
     301                  zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  & 
     302                        &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  ) 
     303                  zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  & 
     304                        &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
     305                  taubot(ji,jj) = 0.5_wp * rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
     306                  ! 
     307            END DO 
     308         END DO 
     309      END IF 
     310 
     311      CALL model%prepare_inputs(real(kt, wp),nyear,nmonth,nday,REAL(nsec_day,wp)) 
     312 
     313      ! Retrieve 3D shortwave and store in etot3 
     314      IF (ln_qsr_spec) THEN 
     315         etot3(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_swr) 
     316      ENDIF 
    293317 
    294318      ! Zero rate array of interface-attached state variables 
     
    296320 
    297321      ! 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) 
     322      DO jj=2,jpjm1 
     323         ! Process bottom (get_bottom_sources increments rather than sets, so zero flux array first) 
    300324         flux = 0._wp 
    301          CALL fabm_do_bottom(model,1,jpi,jj,flux,fabm_st2Da(:,jj,jp_fabm_surface+1:)) 
     325         CALL model%get_bottom_sources(fs_2,fs_jpim1,jj,flux,fabm_st2Da(fs_2:fs_jpim1,jj,jp_fabm_surface+1:)) 
    302326         DO jn=1,jp_fabm 
    303              DO ji=1,jpi 
    304                  ! Divide bottom fluxes by height of bottom layer and add to source terms. 
    305                  ! TODO: is there perhaps an existing variable for fse3t(ji,jj,mbkt(ji,jj))?? 
    306                  tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) = tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) + flux(ji,jn)/fse3t(ji,jj,mbkt(ji,jj)) 
    307              END DO 
    308          END DO 
    309  
    310          ! Process surface (fabm_do_surface increments rather than sets, so zero flux array first) 
     327            ! Divide bottom fluxes by height of bottom layer and add to source terms. 
     328            DO ji=fs_2,fs_jpim1 
     329               tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) = tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) + flux(ji,jn)/fse3t(ji,jj,mbkt(ji,jj)) 
     330            END DO 
     331         END DO 
     332 
     333         ! Process surface (get_surface_sources increments rather than sets, so zero flux array first) 
    311334         flux = 0._wp 
    312          CALL fabm_do_surface(model,1,jpi,jj,flux,fabm_st2Da(:,jj,1:jp_fabm_surface)) 
     335         CALL model%get_surface_sources(fs_2,fs_jpim1,jj,flux,fabm_st2Da(fs_2:fs_jpim1,jj,1:jp_fabm_surface)) 
     336         ! Divide surface fluxes by height of surface layer and add to source terms. 
    313337         DO jn=1,jp_fabm 
    314              ! Divide surface fluxes by height of surface layer and add to source terms. 
    315              tra(:,jj,1,jp_fabm_m1+jn) = tra(:,jj,1,jp_fabm_m1+jn) + flux(:,jn)/fse3t(:,jj,1) 
    316          END DO 
    317       END DO 
    318  
    319       ! Compute interior source terms (NB fabm_do increments rather than sets) 
    320       DO jk=1,jpk 
    321           DO jj=1,jpj 
    322               CALL fabm_do(model,1,jpi,jj,jk,tra(:,jj,jk,jp_fabm0:jp_fabm1)) 
    323           END DO 
    324       END DO 
     338            DO ji=fs_2,fs_jpim1 
     339               tra(ji,jj,1,jp_fabm_m1+jn) = tra(ji,jj,1,jp_fabm_m1+jn) + flux(ji,jn)/fse3t(ji,jj,1) 
     340            END DO 
     341         END DO 
     342      END DO 
     343 
     344      ! Compute interior source terms (NB get_interior_sources increments rather than sets) 
     345      DO jk=1,jpkm1 
     346         DO jj=2,jpjm1 
     347            CALL model%get_interior_sources(fs_2,fs_jpim1,jj,jk,tra(fs_2:fs_jpim1,jj,jk,jp_fabm0:jp_fabm1)) 
     348         END DO 
     349      END DO 
     350 
     351      CALL model%finalize_outputs() 
    325352   END SUBROUTINE compute_fabm 
    326353 
    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 
     354   SUBROUTINE check_state(repair, valid, repaired) 
     355      LOGICAL, INTENT(IN)  :: repair 
     356      LOGICAL, INTENT(OUT) :: valid, repaired 
     357 
     358      INTEGER :: jj, jk 
     359      LOGICAL :: valid_int, valid_sf, valid_bt 
     360 
     361      valid = .TRUE.     ! Whether the model state is valid after this subroutine returns 
     362      repaired = .FALSE. ! Whether the model state has been repaired by this subroutine 
     363      DO jk=1,jpkm1 
     364         DO jj=2,jpjm1 
     365            CALL model%check_interior_state(fs_2, fs_jpim1, jj, jk, repair, valid_int) 
     366            IF (repair .AND. .NOT. valid_int) THEN 
    340367               repair_interior_count = repair_interior_count + 1 
    341                exit_state%repaired = .TRUE. 
     368               repaired = .TRUE. 
    342369            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 
     370            IF (.NOT. (valid_int .OR. repair)) valid = .FALSE. 
     371         END DO 
     372      END DO 
     373      DO jj=2,jpjm1 
     374         CALL model%check_surface_state(fs_2, fs_jpim1, jj, repair, valid_sf) 
     375         IF (repair .AND. .NOT. valid_sf) THEN 
    349376            repair_surface_count = repair_surface_count + 1 
    350             exit_state%repaired = .TRUE. 
     377            repaired = .TRUE. 
    351378         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 
     379         IF (.NOT. (valid_sf .AND. valid_bt) .AND. .NOT. repair) valid = .FALSE. 
     380         CALL model%check_bottom_state(fs_2, fs_jpim1, jj, repair, valid_bt) 
     381         IF (repair .AND. .NOT. valid_bt) THEN 
    355382            repair_bottom_count = repair_bottom_count + 1 
    356             exit_state%repaired = .TRUE. 
     383            repaired = .TRUE. 
    357384         END IF 
    358          IF (.NOT.(valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE. 
    359       END DO 
    360    END FUNCTION 
     385         IF (.NOT. (valid_sf .AND. valid_bt) .AND. .NOT. repair) valid = .FALSE. 
     386      END DO 
     387   END SUBROUTINE 
    361388 
    362389   SUBROUTINE trc_sms_fabm_check_mass() 
    363390      REAL(wp) :: total(SIZE(model%conserved_quantities)) 
    364       INTEGER :: jk,jj,jn 
     391      INTEGER :: ji,jk,jj,jn 
    365392 
    366393      total = 0._wp 
    367394 
    368       DO jk=1,jpk 
    369          DO jj=1,jpj 
    370             CALL fabm_get_conserved_quantities(model,1,jpi,jj,jk,current_total) 
     395      IF (.NOT. started) CALL nemo_fabm_start 
     396 
     397      DO jk=1,jpkm1 
     398         DO jj=2,jpjm1 
     399            CALL model%get_interior_conserved_quantities(fs_2,fs_jpim1,jj,jk,current_total) 
    371400            DO jn=1,SIZE(model%conserved_quantities) 
    372                total(jn) = total(jn) + SUM(cvol(:,jj,jk)*current_total(:,jn)*tmask_i(:,jj)) 
     401               DO ji=fs_2,fs_jpim1 
     402                  total(jn) = total(jn) + cvol(ji,jj,jk) * current_total(ji,jn) * tmask_i(ji,jj) 
     403               END DO 
    373404            END DO 
    374405         END DO 
    375406      END DO 
    376407 
    377       DO jj=1,jpj 
    378          CALL fabm_get_horizontal_conserved_quantities(model,1,jpi,jj,current_total) 
     408      DO jj=2,jpjm1 
     409         CALL model%get_horizontal_conserved_quantities(fs_2,fs_jpim1,jj,current_total) 
    379410         DO jn=1,SIZE(model%conserved_quantities) 
    380             total(jn) = total(jn) + SUM(e1e2t(:,jj)*current_total(:,jn)*tmask_i(:,jj)) 
     411            DO ji=fs_2,fs_jpim1 
     412               total(jn) = total(jn) + e1e2t(ji,jj) * current_total(ji,jn) * tmask_i(ji,jj) 
     413            END DO 
    381414         END DO 
    382415      END DO 
     
    434467 
    435468   INTEGER FUNCTION trc_sms_fabm_alloc() 
    436       INTEGER :: jj,jk,jn 
     469      INTEGER :: jn 
    437470      !!---------------------------------------------------------------------- 
    438471      !!              ***  ROUTINE trc_sms_fabm_alloc  *** 
     
    441474      ! ALLOCATE here the arrays specific to FABM 
    442475      ALLOCATE( lk_rad_fabm(jp_fabm)) 
    443       ALLOCATE( prn(jpi, jpj, jpk)) 
    444       ALLOCATE( rho(jpi, jpj, jpk)) 
    445       ALLOCATE( taubot(jpi, jpj)) 
     476      IF (model%variable_needs_values(fabm_standard_variables%pressure)) ALLOCATE(prn(jpi, jpj, jpk)) 
     477      IF (ALLOCATED(prn) .or. model%variable_needs_values(fabm_standard_variables%density)) ALLOCATE(rho(jpi, jpj, jpk)) 
     478      IF (model%variable_needs_values(fabm_standard_variables%bottom_stress)) ALLOCATE(taubot(jpi, jpj)) 
    446479      ! ALLOCATE( tab(...) , STAT=trc_sms_fabm_alloc ) 
    447480 
     
    452485 
    453486      ! Work array to hold surface and bottom fluxes 
    454       ALLOCATE(flux(jpi,jp_fabm)) 
    455  
    456       ! Work array to hold extinction coefficients 
    457       ALLOCATE(ext(jpi)) 
    458       ext=0._wp 
     487      ALLOCATE(flux(fs_2:fs_jpim1,jp_fabm)) 
    459488 
    460489      ! Allocate work arrays for vertical movement 
    461       ALLOCATE(w_ct(jpi,jpk,jp_fabm)) 
    462       ALLOCATE(w_if(jpk,jp_fabm)) 
    463       ALLOCATE(zwgt_if(jpk,jp_fabm)) 
    464       ALLOCATE(flux_if(jpk,jp_fabm)) 
    465       ALLOCATE(current_total(jpi,SIZE(model%conserved_quantities))) 
     490      ALLOCATE(w_ct(fs_2:fs_jpim1,1:jpkm1,jp_fabm)) 
     491      ALLOCATE(current_total(fs_2:fs_jpim1,SIZE(model%conserved_quantities))) 
    466492#if defined key_trdtrc && defined key_iomput 
    467493      IF( lk_trdtrc ) ALLOCATE(tr_vmv(jpi,jpj,jpk,jp_fabm)) 
     
    474500      ! 
    475501 
    476       ! Make FABM aware of diagnostics that are not needed [not included in output] 
    477       DO jn=1,size(model%diagnostic_variables) 
    478           !model%diagnostic_variables(jn)%save = iom_use(model%diagnostic_variables(jn)%name) 
    479       END DO 
    480       DO jn=1,size(model%horizontal_diagnostic_variables) 
    481           !model%horizontal_diagnostic_variables(jn)%save = iom_use(model%horizontal_diagnostic_variables(jn)%name) 
    482       END DO 
    483  
    484       ! Provide FABM with domain extents [after this, the save attribute of diagnostic variables can no longe change!] 
    485       call fabm_set_domain(model,jpi, jpj, jpk) 
    486  
    487       ! Provide FABM with the vertical indices of the surface and bottom, and the land-sea mask. 
    488       call model%set_bottom_index(mbkt)  ! NB mbkt extents should match dimension lengths provided to fabm_set_domain 
    489       call model%set_surface_index(1) 
    490       call fabm_set_mask(model,tmask,tmask(:,:,1)) ! NB tmask extents should match dimension lengths provided to fabm_set_domain 
    491  
    492       ! Send pointers to state data to FABM 
    493       do jn=1,jp_fabm 
    494          call fabm_link_bulk_state_data(model,jn,trn(:,:,:,jp_fabm_m1+jn)) 
    495       end do 
     502      ! Provide FABM with domain extents 
     503      CALL model%set_domain(jpi, jpj, jpk, rdt) 
     504      CALL model%set_domain_start(fs_2, 2, 1) 
     505      CALL model%set_domain_stop(fs_jpim1, jpjm1, jpkm1) 
     506 
     507      ! Provide FABM with the vertical indices of the bottom, and the land-sea mask. 
     508      CALL model%set_bottom_index(mbkt)  ! NB mbkt extents should match dimension lengths provided to set_domain 
     509      CALL model%set_mask(tmask,tmask(:,:,1)) ! NB tmask extents should match dimension lengths provided to set_domain 
     510 
     511      ! Initialize state and send pointers to state data to FABM 
     512      ! We mask land points in states with zeros, as per with NEMO "convention" 
     513      ! NB we cannot call model%initialize_*_state at this point, because model%start has not been called yet. 
     514      DO jn=1,jp_fabm 
     515         trn(:,:,:,jp_fabm_m1+jn) = model%interior_state_variables(jn)%initial_value * tmask 
     516         CALL model%link_interior_state_data(jn,trn(:,:,:,jp_fabm_m1+jn)) 
     517      END DO 
    496518      DO jn=1,jp_fabm_surface 
    497          CALL fabm_link_surface_state_data(model,jn,fabm_st2Dn(:,:,jn)) 
     519         fabm_st2Dn(:,:,jn) = model%surface_state_variables(jn)%initial_value * tmask(:,:,1) 
     520         CALL model%link_surface_state_data(jn,fabm_st2Dn(:,:,jn)) 
    498521      END DO 
    499522      DO jn=1,jp_fabm_bottom 
    500          CALL fabm_link_bottom_state_data(model,jn,fabm_st2Dn(:,:,jp_fabm_surface+jn)) 
     523         fabm_st2Dn(:,:,jp_fabm_surface+jn) = model%bottom_state_variables(jn)%initial_value * tmask(:,:,1) 
     524         CALL model%link_bottom_state_data(jn,fabm_st2Dn(:,:,jp_fabm_surface+jn)) 
    501525      END DO 
    502526 
    503527      ! Send pointers to environmental data to FABM 
    504       call fabm_link_bulk_data(model,standard_variables%temperature,tsn(:,:,:,jp_tem)) 
    505       call fabm_link_bulk_data(model,standard_variables%practical_salinity,tsn(:,:,:,jp_sal)) 
    506       call fabm_link_bulk_data(model,standard_variables%density,rho(:,:,:)) 
    507       call fabm_link_bulk_data(model,standard_variables%pressure,prn) 
    508       call fabm_link_horizontal_data(model,standard_variables%bottom_stress,taubot(:,:)) 
    509       ! correct target for cell thickness depends on NEMO configuration: 
    510 #ifdef key_vvl 
    511       call fabm_link_bulk_data(model,standard_variables%cell_thickness,e3t_n) 
    512 #else 
    513       call fabm_link_bulk_data(model,standard_variables%cell_thickness,e3t_0) 
    514 #endif 
    515       call fabm_link_horizontal_data(model,standard_variables%latitude,gphit) 
    516       call fabm_link_horizontal_data(model,standard_variables%longitude,glamt) 
    517       call fabm_link_scalar_data(model,standard_variables%number_of_days_since_start_of_the_year,daynumber_in_year) 
    518       call fabm_link_horizontal_data(model,standard_variables%wind_speed,wndm(:,:)) 
    519       call fabm_link_horizontal_data(model,standard_variables%surface_downwelling_shortwave_flux,qsr(:,:)) 
    520       call fabm_link_horizontal_data(model,standard_variables%bottom_depth_below_geoid,bathy(:,:)) 
    521  
    522       swr_id = model%get_bulk_variable_id(standard_variables%downwelling_shortwave_flux) 
     528      CALL model%link_interior_data(fabm_standard_variables%depth, fsdept(:,:,:)) 
     529      CALL model%link_interior_data(fabm_standard_variables%temperature, tsn(:,:,:,jp_tem)) 
     530      CALL model%link_interior_data(fabm_standard_variables%practical_salinity, tsn(:,:,:,jp_sal)) 
     531      IF (ALLOCATED(rho)) CALL model%link_interior_data(fabm_standard_variables%density, rho(:,:,:)) 
     532      IF (ALLOCATED(prn)) CALL model%link_interior_data(fabm_standard_variables%pressure, prn) 
     533      IF (ALLOCATED(taubot)) CALL model%link_horizontal_data(fabm_standard_variables%bottom_stress, taubot(:,:)) 
     534      CALL model%link_interior_data(fabm_standard_variables%cell_thickness, fse3t(:,:,:)) 
     535      CALL model%link_horizontal_data(fabm_standard_variables%latitude, gphit) 
     536      CALL model%link_horizontal_data(fabm_standard_variables%longitude, glamt) 
     537      CALL model%link_scalar(fabm_standard_variables%number_of_days_since_start_of_the_year, daynumber_in_year) 
     538      CALL model%link_horizontal_data(fabm_standard_variables%wind_speed, wndm(:,:)) 
     539      CALL model%link_horizontal_data(fabm_standard_variables%surface_downwelling_shortwave_flux, qsr(:,:)) 
     540      CALL model%link_horizontal_data(fabm_standard_variables%bottom_depth_below_geoid, bathy(:,:)) 
     541      CALL model%link_horizontal_data(fabm_standard_variables%ice_area_fraction, fr_i(:,:)) 
    523542 
    524543      ! Obtain user-specified input variables (read from NetCDF file) 
    525       call link_inputs 
    526       call update_inputs( nit000, .false. ) 
    527  
    528       ! Check whether FABM has all required data 
    529       call fabm_check_ready(model) 
    530  
    531       ! Initialize state 
    532       DO jj=1,jpj 
    533          CALL fabm_initialize_surface_state(model,1,jpi,jj) 
    534          CALL fabm_initialize_bottom_state(model,1,jpi,jj) 
    535       END DO 
    536       DO jk=1,jpk 
    537          DO jj=1,jpj 
    538             CALL fabm_initialize_state(model,1,jpi,jj,jk) 
    539          END DO 
    540       END DO 
     544      CALL link_inputs 
     545      CALL update_inputs(nit000, .FALSE.) 
    541546 
    542547      ! Set mask for negativity corrections to the relevant states 
    543       lk_rad_fabm = .FALSE. 
     548      lk_rad_fabm(:) = .FALSE. 
    544549      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.' 
     550        IF (model%interior_state_variables(jn)%minimum >= 0._wp) THEN 
     551          lk_rad_fabm(jn) = .TRUE. 
     552          IF(lwp) WRITE(numout,*) 'FABM clipping for '//TRIM(model%interior_state_variables(jn)%name)//' activated.' 
    548553        END IF 
    549554      END DO 
    550555 
    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 
    563556 
    564557      ! Copy initial condition for interface-attached state variables to "previous" state field 
     
    566559      fabm_st2Db = fabm_st2Dn 
    567560 
    568       ! Initialise repair counters 
    569       repair_interior_count = 0 
    570       repair_surface_count = 0 
    571       repair_bottom_count = 0 
    572  
    573561   END FUNCTION trc_sms_fabm_alloc 
     562 
     563   SUBROUTINE nemo_fabm_start() 
     564      INTEGER :: jn 
     565 
     566      ! Make FABM aware of diagnostics that are not needed [not included in output] 
     567      ! This works only after iom has completely initialised, because it depends on iom_use 
     568      DO jn=1,size(model%interior_diagnostic_variables) 
     569         model%interior_diagnostic_variables(jn)%save = iom_use(model%interior_diagnostic_variables(jn)%name) & 
     570            .or. iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'_VINT') & 
     571            .or. iom_use(TRIM(model%interior_diagnostic_variables(jn)%name)//'25h') & 
     572            .or. iom_use('top_'//TRIM(model%interior_diagnostic_variables(jn)%name)) & 
     573            .or. iom_use('mid_'//TRIM(model%interior_diagnostic_variables(jn)%name)) & 
     574            .or. iom_use('bot_'//TRIM(model%interior_diagnostic_variables(jn)%name)) 
     575      END DO 
     576      model%interior_diagnostic_variables(jp_fabm_o3ta)%save = .TRUE. 
     577      model%interior_diagnostic_variables(jp_fabm_o3ph)%save = .TRUE. 
     578      model%interior_diagnostic_variables(jp_fabm_o3pc)%save = .TRUE. 
     579      model%interior_diagnostic_variables(jp_fabm_pgrow)%save = .TRUE. 
     580      model%interior_diagnostic_variables(jp_fabm_ploss)%save = .TRUE. 
     581      IF( ln_qsr_spec ) THEN 
     582         model%interior_diagnostic_variables(jp_fabm_swr)%save = .TRUE. 
     583      ENDIF 
     584      IF ( jp_fabm_kd490 /= -1 ) THEN 
     585         model%interior_diagnostic_variables(jp_fabm_kd490)%save = .TRUE. 
     586      ENDIF 
     587      IF ( jp_fabm_xeps /= -1 ) THEN 
     588         model%interior_diagnostic_variables(jp_fabm_xeps)%save = .TRUE. 
     589      ENDIF 
     590      DO jn=1,size(model%horizontal_diagnostic_variables) 
     591         model%horizontal_diagnostic_variables(jn)%save = iom_use(model%horizontal_diagnostic_variables(jn)%name) & 
     592            .or. iom_use(TRIM(model%horizontal_diagnostic_variables(jn)%name)//'25h') 
     593      END DO 
     594 
     595      ! Check whether FABM has all required data 
     596      ! [after this, the save attribute of diagnostic variables can no longer change!] 
     597      CALL model%start() 
     598 
     599      started = .TRUE. 
     600 
     601   END SUBROUTINE 
    574602 
    575603#else 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/trcwri_fabm.F90

    r10270 r13576  
    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/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90

    r12352 r13576  
    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/NEMOGCM/NEMO/TOP_SRC/trcini.F90

    r12432 r13576  
    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/NEMOGCM/NEMO/TOP_SRC/trcnam.F90

    r10162 r13576  
    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.