Changeset 5883


Ignore:
Timestamp:
2015-11-13T08:01:08+01:00 (5 years ago)
Author:
gm
Message:

#1613: vvl by default: TRA/TRC remove optimization associated with linear free surface

Location:
branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM
Files:
1 deleted
41 edited
1 moved

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_cfg

    r5866 r5883  
    257257/ 
    258258!----------------------------------------------------------------------- 
    259 &namhsb       !  Heat and salt budgets 
    260 !----------------------------------------------------------------------- 
     259&namhsb       !  Heat and salt budgets                                  (default F) 
     260!----------------------------------------------------------------------- 
     261   ln_diahsb  = .true.    !  check the heat and salt budgets (T) or not (F) 
    261262/ 
    262263!----------------------------------------------------------------------- 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/CONFIG/SHARED/namelist_ref

    r5866 r5883  
    3333   nn_leapy    =       0   !  Leap year calendar (1) or not (0) 
    3434   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
    35    nn_euler    =       1   !  = 0 : start with forward time step if ln_rstart=T 
    36    nn_rstctl   =       0   !  restart control ==> activated only if ln_rstart=T 
    37                            !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
    38                            !    = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart 
    39                            !    = 2 nn_date0 read in restart  ; nn_it000 : check consistancy between namelist and restart 
    40    cn_ocerst_in  = "restart"   !  suffix of ocean restart name (input) 
    41    cn_ocerst_indir = "."       !  directory from which to read input ocean restarts 
    42    cn_ocerst_out = "restart"   !  suffix of ocean restart name (output) 
    43    cn_ocerst_outdir = "."      !  directory in which to write output ocean restarts 
     35      nn_euler    =    1         !  = 0 : start with forward time step if ln_rstart=T 
     36      nn_rstctl   =    0         !  restart control ==> activated only if ln_rstart=T 
     37                                 !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
     38                                 !    = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart 
     39                                 !    = 2 nn_date0 read in restart  ; nn_it000 : check consistancy between namelist and restart 
     40      cn_ocerst_in  = "restart"  !  suffix of ocean restart name (input) 
     41      cn_ocerst_indir = "."      !  directory from which to read input ocean restarts 
     42      cn_ocerst_out = "restart"  !  suffix of ocean restart name (output) 
     43      cn_ocerst_outdir = "."     !  directory in which to write output ocean restarts 
    4444   nn_istate   =       0   !  output the initial state (1) or not (0) 
    4545   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
     
    5050   ln_mskland  = .false.   !  mask land points in NetCDF outputs (costly: + ~15%) 
    5151   ln_cfmeta   = .false.   !  output additional data to netCDF files required for compliance with the CF metadata standard 
    52    ln_clobber  = .true.   !  clobber (overwrite) an existing file 
     52   ln_clobber  = .true.    !  clobber (overwrite) an existing file 
    5353   nn_chunksz  =       0   !  chunksize (bytes) for NetCDF file (works only with iom_nf90 routines) 
    5454/ 
     
    6565! 
    6666!----------------------------------------------------------------------- 
    67 &namcfg     !   parameters of the configuration 
    68 !----------------------------------------------------------------------- 
    69    cp_cfg      =  "default"            !  name of the configuration 
    70    cp_cfz      =  "no zoom"            !  name of the zoom of configuration 
    71    jp_cfg      =       0               !  resolution of the configuration 
    72    jpidta      =      10               !  1st lateral dimension ( >= jpi ) 
    73    jpjdta      =      12               !  2nd    "         "    ( >= jpj ) 
    74    jpkdta      =      31               !  number of levels      ( >= jpk ) 
    75    jpiglo      =      10               !  1st dimension of global domain --> i =jpidta 
    76    jpjglo      =      12               !  2nd    -                  -    --> j =jpjdta 
    77    jpizoom     =       1               !  left bottom (i,j) indices of the zoom 
    78    jpjzoom     =       1               !  in data domain indices 
    79    jperio      =       0               !  lateral cond. type (between 0 and 6) 
    80                                        !  = 0 closed                 ;   = 1 cyclic East-West 
    81                                        !  = 2 equatorial symmetric   ;   = 3 North fold T-point pivot 
    82                                        !  = 4 cyclic East-West AND North fold T-point pivot 
    83                                        !  = 5 North fold F-point pivot 
    84                                        !  = 6 cyclic East-West AND North fold F-point pivot 
    85    ln_use_jattr = .false.              !  use (T) the file attribute: open_ocean_jstart, if present 
    86                                        !  in netcdf input files, as the start j-row for reading 
     67&namcfg        !   parameters of the configuration 
     68!----------------------------------------------------------------------- 
     69   cp_cfg      = "default" !  name of the configuration 
     70   cp_cfz      = "no zoom" !  name of the zoom of configuration 
     71   jp_cfg      =      0    !  resolution of the configuration 
     72   jpidta      =     10    !  1st lateral dimension ( >= jpi ) 
     73   jpjdta      =     12    !  2nd    "         "    ( >= jpj ) 
     74   jpkdta      =     31    !  number of levels      ( >= jpk ) 
     75   jpiglo      =     10    !  1st dimension of global domain --> i =jpidta 
     76   jpjglo      =     12    !  2nd    -                  -    --> j =jpjdta 
     77   jpizoom     =      1    !  left bottom (i,j) indices of the zoom 
     78   jpjzoom     =      1    !  in data domain indices 
     79   jperio      =      0    !  lateral cond. type (between 0 and 6) 
     80                                 !  = 0 closed                 ;   = 1 cyclic East-West 
     81                                 !  = 2 equatorial symmetric   ;   = 3 North fold T-point pivot 
     82                                 !  = 4 cyclic East-West AND North fold T-point pivot 
     83                                 !  = 5 North fold F-point pivot 
     84                                 !  = 6 cyclic East-West AND North fold F-point pivot 
     85   ln_use_jattr = .false.  !  use (T) the file attribute: open_ocean_jstart, if present 
     86                           !  in netcdf input files, as the start j-row for reading 
    8787/ 
    8888!----------------------------------------------------------------------- 
     
    236236!----------------------------------------------------------------------- 
    237237   nn_fsbc     = 5         !  frequency of surface boundary condition computation 
    238                            !     (also = the frequency of sea-ice model call) 
     238                           !     (also = the frequency of sea-ice & iceberg model call) 
     239                     ! Type of air-sea fluxes  
    239240   ln_ana      = .false.   !  analytical formulation                    (T => fill namsbc_ana ) 
    240241   ln_flx      = .false.   !  flux formulation                          (T => fill namsbc_flx ) 
     
    242243   ln_blk_core = .true.    !  CORE bulk formulation                     (T => fill namsbc_core) 
    243244   ln_blk_mfs  = .false.   !  MFS bulk formulation                      (T => fill namsbc_mfs ) 
     245                     ! Type of coupling (Ocean/Ice/Atmosphere) : 
    244246   ln_cpl      = .false.   !  atmosphere coupled   formulation          ( requires key_oasis3 ) 
    245247   ln_mixcpl   = .false.   !  forced-coupled mixed formulation          ( requires key_oasis3 ) 
     
    248250                           !  =1 opa-sas OASIS coupling: multi executable configuration, OPA component 
    249251                           !  =2 opa-sas OASIS coupling: multi executable configuration, SAS component  
    250    ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr ) 
     252   nn_limflx = -1          !  LIM3 Multi-category heat flux formulation (use -1 if LIM3 is not used) 
     253                           !  =-1  Use per-category fluxes, bypass redistributor, forced mode only, not yet implemented coupled 
     254                           !  = 0  Average per-category fluxes (forced and coupled mode) 
     255                           !  = 1  Average and redistribute per-category fluxes, forced mode only, not yet implemented coupled 
     256                           !  = 2  Redistribute a single flux over categories (coupled mode only) 
     257                     ! Sea-ice : 
    251258   nn_ice      = 2         !  =0 no ice boundary condition   , 
    252259                           !  =1 use observed ice-cover      , 
     
    255262                           !  =1 levitating ice with mass and salt exchange but no presure effect 
    256263                           !  =2 embedded sea-ice (full salt and mass exchanges and pressure) 
    257    ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    258    ln_rnf      = .true.    !  runoffs                                   (T   => fill namsbc_rnf) 
     264                     ! Misc. options of sbc :  
     265   ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr ) 
     266      ln_dm2dc = .false.      !  daily mean to diurnal cycle on short wave 
     267   ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf) 
     268   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr) 
     269   nn_fwb      = 2         !  FreshWater Budget: =0 unchecked 
     270                           !     =1 global mean of e-p-r set to zero at each time step 
     271                           !     =2 annual global mean of e-p-r set to zero 
     272   ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr ) 
    259273   nn_isf      = 0         !  ice shelf melting/freezing                (/=0 => fill namsbc_isf) 
    260274                           !  0 =no isf                  1 = presence of ISF 
     
    262276                           !  4 = ISF fwf specified 
    263277                           !  option 1 and 4 need ln_isfcav = .true. (domzgr) 
    264    ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr) 
    265    nn_fwb      = 2         !  FreshWater Budget: =0 unchecked 
    266                            !     =1 global mean of e-p-r set to zero at each time step 
    267                            !     =2 annual global mean of e-p-r set to zero 
    268    ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave) 
    269    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave) 
    270    ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave) 
     278   ln_wave = .false.       !  coupling with surface wave                    (T => fill namsbc_wave) 
    271279   nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) , 
    272280                           !  =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field) 
    273    nn_limflx = -1          !  LIM3 Multi-category heat flux formulation (use -1 if LIM3 is not used) 
    274                            !  =-1  Use per-category fluxes, bypass redistributor, forced mode only, not yet implemented coupled 
    275                            !  = 0  Average per-category fluxes (forced and coupled mode) 
    276                            !  = 1  Average and redistribute per-category fluxes, forced mode only, not yet implemented coupled 
    277                            !  = 2  Redistribute a single flux over categories (coupled mode only) 
    278281/ 
    279282!----------------------------------------------------------------------- 
     
    406409 
    407410   cn_dir      = './'      !  root directory for the location of the runoff files 
    408    ln_traqsr   = .true.    !  Light penetration (T) or not (F) 
    409411   ln_qsr_rgb  = .true.    !  RGB (Red-Green-Blue) light penetration 
    410412   ln_qsr_2bd  = .false.   !  2 bands              light penetration 
     
    11521154/ 
    11531155!----------------------------------------------------------------------- 
    1154 &namhsb       !  Heat and salt budgets 
     1156&namhsb       !  Heat and salt budgets                                  (default F) 
    11551157!----------------------------------------------------------------------- 
    11561158   ln_diahsb  = .false.    !  check the heat and salt budgets (T) or not (F) 
     
    12671269   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    12681270! 
    1269    cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files 
    1270 / 
     1271   cn_dir_cdg  = './'      !  root directory for the location of drag coefficient files 
     1272   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model 
     1273   ln_sdw  = .false.       !  Computation of 3D stokes drift                
     1274/ 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r5866 r5883  
    869869            ! 
    870870            sshb(:,:) = sshn(:,:)                        ! Update before fields 
    871             ! 
    872             IF( .NOT.ln_linssh ) THEN 
    873                DO jk = 1, jpk 
    874                   e3t_b(:,:,jk) = e3t_n(:,:,jk) 
    875                END DO 
    876             ENDIF 
     871            e3t_b(:,:,:) = e3t_n(:,:,:) 
     872!!gm why not e3u_b, e3v_b, gdept_b ???? 
    877873            ! 
    878874            DEALLOCATE( ssh_bkg    ) 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r5866 r5883  
    113113            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
    114114         END IF 
     115!!gm 
     116!!gm   riceload should be added in both ln_linssh=T or F, no? 
     117!!gm 
    115118      END IF 
    116119      !                                          
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5866 r5883  
    103103      !              !==  time varying part of coordinate system  ==! 
    104104      ! 
    105       IF( .NOT.ln_linssh ) THEN                 ! time varying : initialize before/now/after variables 
    106          ! 
    107          CALL dom_vvl_init  
    108          ! 
    109       ELSE                                      ! Fix in time : set to the reference one for all 
     105      IF( ln_linssh ) THEN          ! Fix in time : set to the reference one for all 
    110106         !       before        !          now          !       after         ! 
    111107         ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points 
     
    134130         ! 
    135131         CALL wrk_dealloc( jpi,jpj,   z1_hu_0, z1_hv_0 ) 
     132         ! 
     133      ELSE                         ! time varying : initialize before/now/after variables 
     134         ! 
     135         CALL dom_vvl_init  
     136         ! 
    136137      ENDIF 
    137138      ! 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5866 r5883  
    267267      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling. 
    268268      !!---------------------------------------------------------------------- 
    269       REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 
    270       REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv 
    271       !! * Arguments 
    272       INTEGER, INTENT( in )                  :: kt                    ! time step 
    273       INTEGER, INTENT( in ), OPTIONAL        :: kcall                 ! optional argument indicating call sequence 
    274       !! * Local declarations 
     269      INTEGER, INTENT( in )           ::   kt      ! time step 
     270      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence 
     271      ! 
    275272      INTEGER                                :: ji, jj, jk            ! dummy loop indices 
    276273      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers 
     
    278275      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars 
    279276      LOGICAL                                :: ll_do_bclinic         ! temporary logical 
    280       !!---------------------------------------------------------------------- 
     277      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 
     278      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv 
     279      !!---------------------------------------------------------------------- 
     280      ! 
     281      IF( ln_linssh )   RETURN      ! No calculation in linear free surface 
    281282      ! 
    282283      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt') 
     
    579580      REAL(wp) ::   zcoef        ! local scalar 
    580581      !!---------------------------------------------------------------------- 
    581  
     582      ! 
     583      IF( ln_linssh )   RETURN      ! No calculation in linear free surface 
     584      ! 
    582585      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp') 
    583586      ! 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5866 r5883  
    122122            ENDIF 
    123123         ENDIF 
    124          !    
     124         ! 
     125!!gm This is to be changed !!!! 
    125126         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here 
    126127         IF( .NOT.ln_linssh ) THEN 
     
    129130            END DO 
    130131         ENDIF 
     132!!gm  
    131133         !  
    132134      ENDIF 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90

    r5866 r5883  
    1010 
    1111   !!---------------------------------------------------------------------- 
    12    !!   dyn_adv_cen2       : flux form momentum advection (ln_dynadv_cen2=T) 
    13    !!                        trends using a 2nd order centred scheme   
     12   !!   dyn_adv_cen2  : flux form momentum advection (ln_dynadv_cen2=T) using a 2nd order centred scheme   
    1413   !!---------------------------------------------------------------------- 
    1514   USE oce            ! ocean dynamics and tracers 
     
    6766      ENDIF 
    6867      ! 
    69       IF( l_trddyn ) THEN           ! Save ua and va trends 
     68      IF( l_trddyn ) THEN           ! trends: store the input trends 
    7069         zfu_uw(:,:,:) = ua(:,:,:) 
    7170         zfv_vw(:,:,:) = va(:,:,:) 
    7271      ENDIF 
    73  
    74       !                                      ! ====================== ! 
    75       !                                      !  Horizontal advection  ! 
    76       DO jk = 1, jpkm1                       ! ====================== ! 
    77          !                                         ! horizontal volume fluxes 
     72      ! 
     73      !                             !==  Horizontal advection  ==! 
     74      ! 
     75      DO jk = 1, jpkm1                    ! horizontal transport 
    7876         zfu(:,:,jk) = 0.25 * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    7977         zfv(:,:,jk) = 0.25 * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
    80          ! 
    81          DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point 
     78         DO jj = 1, jpjm1                 ! horizontal momentum fluxes (at T- and F-point) 
    8279            DO ji = 1, fs_jpim1   ! vector opt. 
    83                zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj  ,jk) ) 
    84                zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) ) 
    85                zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) ) 
    86                zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) ) 
     80               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj  ,jk) ) 
     81               zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) ) 
     82               zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) ) 
     83               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) ) 
    8784            END DO 
    8885         END DO 
    89          DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes 
     86         DO jj = 2, jpjm1                 ! divergence of horizontal momentum fluxes 
    9087            DO ji = fs_2, fs_jpim1   ! vector opt. 
    9188               zbu = e1e2u(ji,jj) * e3u_n(ji,jj,jk) 
    9289               zbv = e1e2v(ji,jj) * e3v_n(ji,jj,jk) 
    9390               ! 
    94                ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    & 
    95                   &                           + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu 
    96                va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)    & 
    97                   &                           + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv 
     91               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
     92                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) / zbu 
     93               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
     94                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) / zbv 
    9895            END DO 
    9996         END DO 
    10097      END DO 
    10198      ! 
    102       IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic 
     99      IF( l_trddyn ) THEN           ! trends: send trend to trddyn for diagnostic 
    103100         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 
    104101         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 
     
    108105      ENDIF 
    109106      ! 
    110  
    111       !                                      ! ==================== ! 
    112       !                                      !  Vertical advection  ! 
    113       DO jk = 1, jpkm1                       ! ==================== ! 
    114          !                                         ! Vertical volume fluxesÊ 
    115          zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk) 
    116          ! 
    117          IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                    
    118             zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero 
    119             zfv_vw(:,:,jpk) = 0.e0 
    120             !                                           ! Surface value : 
    121             IF( .NOT.ln_linssh ) THEN                        ! variable volume : flux set to zero 
    122                zfu_uw(:,:, 1 ) = 0._wp     
    123                zfv_vw(:,:, 1 ) = 0._wp 
    124             ELSE                                             ! constant volume : advection through the surface 
    125                DO jj = 2, jpjm1 
    126                   DO ji = fs_2, fs_jpim1 
    127                      zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1) 
    128                      zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1) 
    129                   END DO 
    130                END DO 
    131             ENDIF 
    132          ELSE                                      ! interior fluxes 
    133             DO jj = 2, jpjm1 
    134                DO ji = fs_2, fs_jpim1   ! vector opt. 
    135                   zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
    136                   zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
    137                END DO 
     107      !                             !==  Vertical advection  ==! 
     108      ! 
     109      DO jj = 2, jpjm1                    ! surface/bottom advective fluxes set to zero 
     110         DO ji = fs_2, fs_jpim1 
     111            zfu_uw(ji,jj,jpk) = 0._wp   ;   zfv_vw(jj,jj,jpk) = 0._wp 
     112            zfu_uw(ji,jj, 1 ) = 0._wp   ;   zfv_vw(jj,jj, 1 ) = 0._wp 
     113         END DO 
     114      END DO 
     115      IF( ln_linssh ) THEN                ! linear free surface: advection through the surface 
     116         DO jj = 2, jpjm1 
     117            DO ji = fs_2, fs_jpim1 
     118               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1) 
     119               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1) 
    138120            END DO 
    139          ENDIF 
    140       END DO 
    141       DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence 
    142          DO jj = 2, jpjm1  
     121         END DO 
     122      ENDIF 
     123      DO jk = 2, jpkm1                    ! interior advective fluxes 
     124         DO jj = 2, jpjm1                       ! 1/4 * Vertical transport 
     125            DO ji = fs_2, fs_jpim1 
     126               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     127            END DO 
     128         END DO 
     129         DO jj = 2, jpjm1 
    143130            DO ji = fs_2, fs_jpim1   ! vector opt. 
    144                ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    & 
    145                   &  / ( e1e2u(ji,jj) * e3u_n(ji,jj,jk) ) 
    146                va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    & 
    147                   &  / ( e1e2v(ji,jj) * e3v_n(ji,jj,jk) ) 
     131               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
     132               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
    148133            END DO 
    149134         END DO 
    150135      END DO 
    151136      ! 
    152       IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic 
     137      DO jk = 1, jpkm1                    ! divergence of vertical momentum flux divergence 
     138         DO jj = 2, jpjm1  
     139            DO ji = fs_2, fs_jpim1   ! vector opt. 
     140               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
     141               va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     142            END DO 
     143         END DO 
     144      END DO 
     145      ! 
     146      IF( l_trddyn ) THEN                 ! trends: send trend to trddyn for diagnostic 
    153147         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 
    154148         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 
    155149         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 
    156150      ENDIF 
    157       !                                            ! Control print 
     151      !                                   ! Control print 
    158152      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' cen2 adv - Ua: ', mask1=umask,   & 
    159153         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' ) 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90

    r5866 r5883  
    100100      zlu_uv(:,:,:,:) = 0._wp  
    101101      zlv_vu(:,:,:,:) = 0._wp  
    102  
    103       IF( l_trddyn ) THEN           ! Save ua and va trends 
     102      ! 
     103      IF( l_trddyn ) THEN           ! trends: store the input trends 
    104104         zfu_uw(:,:,:) = ua(:,:,:) 
    105105         zfv_vw(:,:,:) = va(:,:,:) 
    106106      ENDIF 
    107  
    108107      !                                      ! =========================== ! 
    109108      DO jk = 1, jpkm1                       !  Laplacian of the velocity  ! 
     
    115114         DO jj = 2, jpjm1                          ! laplacian 
    116115            DO ji = fs_2, fs_jpim1   ! vector opt. 
    117                ! 
    118116               zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj  ,jk) - 2.*ub (ji,jj,jk) + ub (ji-1,jj  ,jk) ) * umask(ji,jj,jk) 
    119117               zlv_vv(ji,jj,jk,1) = ( vb (ji  ,jj+1,jk) - 2.*vb (ji,jj,jk) + vb (ji  ,jj-1,jk) ) * vmask(ji,jj,jk) 
     
    136134      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. ) 
    137135      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. )  
    138        
     136      ! 
    139137      !                                      ! ====================== ! 
    140138      !                                      !  Horizontal advection  ! 
     
    149147               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) ) 
    150148               ! 
    151                IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1) 
    152                ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1) 
    153                ENDIF 
    154                IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1) 
    155                ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1) 
     149               IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1) 
     150               ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1) 
     151               ENDIF 
     152               IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1) 
     153               ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1) 
    156154               ENDIF 
    157155               ! 
     
    165163               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) ) 
    166164               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) ) 
    167                IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1) 
    168                ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1) 
    169                ENDIF 
    170                IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1) 
    171                ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1) 
     165               IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1) 
     166               ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1) 
     167               ENDIF 
     168               IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1) 
     169               ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1) 
    172170               ENDIF 
    173171               ! 
     
    190188         END DO 
    191189      END DO 
    192       IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic 
     190      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic 
    193191         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 
    194192         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 
     
    197195         zfv_t(:,:,:) = va(:,:,:) 
    198196      ENDIF 
    199  
    200197      !                                      ! ==================== ! 
    201198      !                                      !  Vertical advection  ! 
    202       DO jk = 1, jpkm1                       ! ==================== ! 
    203          !                                         ! Vertical volume fluxesÊ 
    204          zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk) 
    205          ! 
    206          IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                    
    207             zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero 
    208             zfv_vw(:,:,jpk) = 0.e0 
    209             !                                           ! Surface value : 
    210             IF( .NOT.ln_linssh ) THEN                        ! variable volume : flux set to zero 
    211                zfu_uw(:,:, 1 ) = 0._wp 
    212                zfv_vw(:,:, 1 ) = 0._wp 
    213             ELSE                                             ! constant volume : advection through the surface 
    214                DO jj = 2, jpjm1 
    215                   DO ji = fs_2, fs_jpim1 
    216                      zfu_uw(ji,jj, 1 ) = 2._wp * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1) 
    217                      zfv_vw(ji,jj, 1 ) = 2._wp * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1) 
    218                   END DO 
    219                END DO 
    220             ENDIF 
    221          ELSE                                      ! interior fluxes 
    222             DO jj = 2, jpjm1 
    223                DO ji = fs_2, fs_jpim1   ! vector opt. 
    224                   zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
    225                   zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
    226                END DO 
    227             END DO 
    228          ENDIF 
    229       END DO 
    230       DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence 
    231          DO jj = 2, jpjm1  
    232             DO ji = fs_2, fs_jpim1   ! vector opt. 
    233                ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    & 
    234                   &  / ( e1e2u(ji,jj) * e3u_n(ji,jj,jk) ) 
    235                va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    & 
    236                   &  / ( e1e2v(ji,jj) * e3v_n(ji,jj,jk) ) 
    237             END DO 
    238          END DO 
    239       END DO 
    240       ! 
    241       IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic 
     199      !                                      ! ==================== ! 
     200      DO jj = 2, jpjm1                             ! surface/bottom advective fluxes set to zero                   
     201         DO ji = fs_2, fs_jpim1 
     202            zfu_uw(ji,jj,jpk) = 0._wp 
     203            zfv_vw(ji,jj,jpk) = 0._wp 
     204            zfu_uw(ji,jj, 1 ) = 0._wp 
     205            zfv_vw(ji,jj, 1 ) = 0._wp 
     206         END DO 
     207      END DO 
     208      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface 
     209         DO jj = 2, jpjm1 
     210            DO ji = fs_2, fs_jpim1 
     211               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1) 
     212               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1) 
     213            END DO 
     214         END DO 
     215      ENDIF 
     216      DO jk = 2, jpkm1                          ! interior fluxes 
     217         DO jj = 2, jpjm1 
     218            DO ji = fs_2, fs_jpim1 
     219               zfw(ji,jj,jk) = 0.25 * e1e2t(ji,jj) * wn(ji,jj,jk) 
     220            END DO 
     221         END DO 
     222         DO jj = 2, jpjm1 
     223            DO ji = fs_2, fs_jpim1   ! vector opt. 
     224               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
     225               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
     226            END DO 
     227         END DO 
     228      END DO 
     229      DO jk = 1, jpkm1                          ! divergence of vertical momentum flux divergence 
     230         DO jj = 2, jpjm1 
     231            DO ji = fs_2, fs_jpim1   ! vector opt. 
     232               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
     233               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     234            END DO 
     235         END DO 
     236      END DO 
     237      ! 
     238      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic 
    242239         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 
    243240         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 
    244241         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 
    245242      ENDIF 
    246       !                                            ! Control print 
     243      !                                         ! Control print 
    247244      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   & 
    248245         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' ) 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

    r5845 r5883  
    6363        zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
    6464 
    65         IF( l_trddyn )   THEN                      ! temporary save of ua and va trends 
    66            CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
     65        IF( l_trddyn ) THEN      ! trends: store the input trends 
     66           CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv ) 
    6767           ztrdu(:,:,:) = ua(:,:,:) 
    6868           ztrdv(:,:,:) = va(:,:,:) 
     
    8080           END DO 
    8181        END DO 
    82          
    83         IF ( ln_isfcav ) THEN 
     82        ! 
     83        IF( ln_isfcav ) THEN        ! ocean cavities 
    8484           DO jj = 2, jpjm1 
    8585              DO ji = 2, jpim1 
     
    9797           END DO 
    9898        END IF 
    99  
    10099        ! 
    101         IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
     100        IF( l_trddyn ) THEN      ! trends: send trends to trddyn for further diagnostics 
    102101           ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    103102           ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    104103           CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 
    105            CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
     104           CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv ) 
    106105        ENDIF 
    107106        !                                          ! print mean trends (used for debugging) 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap_blp.F90

    r5866 r5883  
    44   !! Ocean dynamics:  lateral viscosity trend (laplacian and bilaplacian) 
    55   !!====================================================================== 
    6    !! History :  OPA  ! 1990-09 (G. Madec) Original code 
    7    !!            4.0  ! 1991-11 (G. Madec) 
    8    !!            6.0  ! 1996-01 (G. Madec) statement function for e3 and ahm 
    9    !!   NEMO     1.0  ! 2002-06 (G. Madec)  F90: Free form and module 
    10    !!             -   ! 2004-08 (C. Talandier) New trends organization 
    11    !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification, 
    12    !!                 !                                  add velocity dependent coefficient and optional read in file 
     6   !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian 
    137   !!---------------------------------------------------------------------- 
    148 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90

    r5429 r5883  
    44   !! Ocean        : lateral boundary conditions 
    55   !!===================================================================== 
    6    !! History :  OPA  ! 1997-06  (G. Madec)     Original code 
    7    !!   NEMO     1.0  ! 2002-09  (G. Madec)     F90: Free form and module 
     6   !! History :  OPA  ! 1997-06  (G. Madec)  Original code 
     7   !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module 
    88   !!            3.2  ! 2009-03  (R. Benshila)  External north fold treatment   
    9    !!            3.5  ! 2012     (S.Mocavero, I. Epicoco) Add 'lbc_bdy_lnk'  
    10    !!                            and lbc_obc_lnk' routine to optimize   
    11    !!                            the BDY/OBC communications 
    12    !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  add a C1D case   
     9   !!            3.5  ! 2012     (S.Mocavero, I. Epicoco) optimization of BDY comm. via lbc_bdy_lnk and lbc_obc_lnk 
     10   !!            3.4  ! 2012-12  (R. Bourdalle-Badie, G. Reffray)  add a C1D case   
    1311   !!---------------------------------------------------------------------- 
    1412#if defined key_mpp_mpi 
     
    2220   USE lib_mpp          ! distributed memory computing library 
    2321 
    24  
    2522   INTERFACE lbc_lnk_multi 
    2623      MODULE PROCEDURE mpp_lnk_2d_9 
    2724   END INTERFACE 
    28  
     25   ! 
    2926   INTERFACE lbc_lnk 
    3027      MODULE PROCEDURE mpp_lnk_3d_gather, mpp_lnk_3d, mpp_lnk_2d 
    3128   END INTERFACE 
    32  
     29   ! 
    3330   INTERFACE lbc_bdy_lnk 
    3431      MODULE PROCEDURE mpp_lnk_bdy_2d, mpp_lnk_bdy_3d 
    3532   END INTERFACE 
    36  
     33   ! 
    3734   INTERFACE lbc_lnk_e 
    3835      MODULE PROCEDURE mpp_lnk_2d_e 
    3936   END INTERFACE 
    40  
     37   ! 
    4138   INTERFACE lbc_lnk_icb 
    4239      MODULE PROCEDURE mpp_lnk_2d_icb 
    4340   END INTERFACE 
    4441 
    45    PUBLIC lbc_lnk       ! ocean lateral boundary conditions 
    46    PUBLIC lbc_lnk_multi ! modified ocean lateral boundary conditions 
    47    PUBLIC lbc_lnk_e 
    48    PUBLIC lbc_bdy_lnk   ! ocean lateral BDY boundary conditions 
    49    PUBLIC lbc_lnk_icb 
     42   PUBLIC   lbc_lnk       ! ocean lateral boundary conditions 
     43   PUBLIC   lbc_lnk_multi ! modified ocean lateral boundary conditions 
     44   PUBLIC   lbc_lnk_e     ! 
     45   PUBLIC   lbc_bdy_lnk   ! ocean lateral BDY boundary conditions 
     46   PUBLIC   lbc_lnk_icb   ! 
    5047 
    5148   !!---------------------------------------------------------------------- 
     
    5451   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    5552   !!---------------------------------------------------------------------- 
    56  
    5753#else 
    5854   !!---------------------------------------------------------------------- 
     
    7571      MODULE PROCEDURE lbc_lnk_3d_gather, lbc_lnk_3d, lbc_lnk_2d 
    7672   END INTERFACE 
    77  
     73   ! 
    7874   INTERFACE lbc_lnk_e 
    7975      MODULE PROCEDURE lbc_lnk_2d_e 
    8076   END INTERFACE 
    81  
     77   ! 
    8278   INTERFACE lbc_bdy_lnk 
    8379      MODULE PROCEDURE lbc_bdy_lnk_2d, lbc_bdy_lnk_3d 
    8480   END INTERFACE 
    85  
     81   ! 
    8682   INTERFACE lbc_lnk_icb 
    8783      MODULE PROCEDURE lbc_lnk_2d_e 
     
    8985 
    9086   PUBLIC   lbc_lnk       ! ocean/ice  lateral boundary conditions 
    91    PUBLIC   lbc_lnk_e  
     87   PUBLIC   lbc_lnk_e     ! 
    9288   PUBLIC   lbc_bdy_lnk   ! ocean lateral BDY boundary conditions 
    93    PUBLIC   lbc_lnk_icb 
     89   PUBLIC   lbc_lnk_icb   ! 
    9490    
    9591   !!---------------------------------------------------------------------- 
     
    230226         ! this is in mpp case. In this module, just do nothing 
    231227      ELSE 
    232          ! 
    233228         !                                     !  East-West boundaries 
    234229         !                                     ! ====================== 
     
    249244            ! 
    250245         END SELECT 
    251          ! 
    252246         !                                     ! North-South boundaries 
    253247         !                                     ! ====================== 
     
    287281   END SUBROUTINE lbc_lnk_3d 
    288282 
     283 
    289284   SUBROUTINE lbc_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval ) 
    290285      !!--------------------------------------------------------------------- 
     
    316311         ! this is in mpp case. In this module, just do nothing 
    317312      ELSE       
    318          ! 
    319313         !                                     ! East-West boundaries 
    320314         !                                     ! ==================== 
     
    335329            ! 
    336330         END SELECT 
    337          ! 
    338331         !                                     ! North-South boundaries 
    339332         !                                     ! ====================== 
     
    375368#endif 
    376369 
    377  
    378370   SUBROUTINE lbc_bdy_lnk_3d( pt3d, cd_type, psgn, ib_bdy ) 
    379371      !!--------------------------------------------------------------------- 
     
    381373      !! 
    382374      !! ** Purpose :   wrapper rountine to 'lbc_lnk_3d'. This wrapper is used 
    383       !!                to maintain the same interface with regards to the mpp 
    384       !case 
    385       !! 
    386       !!---------------------------------------------------------------------- 
    387       CHARACTER(len=1)                , INTENT(in   )           ::   cd_type   ! nature of pt3d grid-points 
    388       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::   pt3d      ! 3D array on which the lbc is applied 
    389       REAL(wp)                        , INTENT(in   )           ::   psgn      ! control of the sign  
    390       INTEGER                                                   ::   ib_bdy    ! BDY boundary set 
    391       !! 
     375      !!              to maintain the same interface with regards to the mpp case 
     376      !! 
     377      !!---------------------------------------------------------------------- 
     378      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type   ! nature of pt3d grid-points 
     379      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pt3d      ! 3D array on which the lbc is applied 
     380      REAL(wp)                        , INTENT(in   ) ::   psgn      ! control of the sign  
     381      INTEGER                         , INTENT(in   ) ::   ib_bdy    ! BDY boundary set 
     382      !!---------------------------------------------------------------------- 
     383      ! 
    392384      CALL lbc_lnk_3d( pt3d, cd_type, psgn) 
    393  
     385      ! 
    394386   END SUBROUTINE lbc_bdy_lnk_3d 
    395387 
     388 
    396389   SUBROUTINE lbc_bdy_lnk_2d( pt2d, cd_type, psgn, ib_bdy ) 
    397390      !!--------------------------------------------------------------------- 
     
    399392      !! 
    400393      !! ** Purpose :   wrapper rountine to 'lbc_lnk_3d'. This wrapper is used 
    401       !!                to maintain the same interface with regards to the mpp 
    402       !case 
    403       !! 
    404       !!---------------------------------------------------------------------- 
    405       CHARACTER(len=1)                , INTENT(in   )           ::   cd_type   ! nature of pt3d grid-points 
    406       REAL(wp), DIMENSION(jpi,jpj),     INTENT(inout)           ::   pt2d      ! 3D array on which the lbc is applied 
    407       REAL(wp)                        , INTENT(in   )           ::   psgn      ! control of the sign  
    408       INTEGER                                                   ::   ib_bdy    ! BDY boundary set 
    409       !! 
     394      !!              to maintain the same interface with regards to the mpp case 
     395      !! 
     396      !!---------------------------------------------------------------------- 
     397      CHARACTER(len=1)            , INTENT(in   ) ::   cd_type   ! nature of pt3d grid-points 
     398      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pt2d      ! 3D array on which the lbc is applied 
     399      REAL(wp)                    , INTENT(in   ) ::   psgn      ! control of the sign  
     400      INTEGER                     , INTENT(in   ) ::   ib_bdy    ! BDY boundary set 
     401      !!---------------------------------------------------------------------- 
     402      ! 
    410403      CALL lbc_lnk_2d( pt2d, cd_type, psgn) 
    411  
     404      ! 
    412405   END SUBROUTINE lbc_bdy_lnk_2d 
    413406 
     
    426419      !!                             for closed boundaries. 
    427420      !!---------------------------------------------------------------------- 
    428       CHARACTER(len=1)            , INTENT(in   )           ::   cd_type   ! nature of pt3d grid-points 
    429       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pt2d      ! 2D array on which the lbc is applied 
    430       REAL(wp)                    , INTENT(in   )           ::   psgn      ! control of the sign  
    431       INTEGER                     , INTENT(in   )           ::   jpri      ! size of extra halo (not needed in non-mpp) 
    432       INTEGER                     , INTENT(in   )           ::   jprj      ! size of extra halo (not needed in non-mpp) 
    433       !!---------------------------------------------------------------------- 
    434  
     421      CHARACTER(len=1)            , INTENT(in   ) ::   cd_type   ! nature of pt3d grid-points 
     422      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pt2d      ! 2D array on which the lbc is applied 
     423      REAL(wp)                    , INTENT(in   ) ::   psgn      ! control of the sign  
     424      INTEGER                     , INTENT(in   ) ::   jpri      ! size of extra halo (not needed in non-mpp) 
     425      INTEGER                     , INTENT(in   ) ::   jprj      ! size of extra halo (not needed in non-mpp) 
     426      !!---------------------------------------------------------------------- 
     427      ! 
    435428      CALL lbc_lnk_2d( pt2d, cd_type, psgn ) 
    436429      !     
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90

    r5836 r5883  
    2323   !!                          the mppobc routine to optimize the BDY and OBC communications 
    2424   !!            3.5  !  2013  ( C. Ethe, G. Madec ) message passing arrays as local variables  
    25    !!            3.5  !  2013 (S.Mocavero, I.Epicoco - CMCC) north fold optimizations 
     25   !!            3.5  !  2013  (S.Mocavero, I.Epicoco - CMCC) north fold optimizations 
    2626   !!---------------------------------------------------------------------- 
    2727 
     
    26622662   END SUBROUTINE mpp_lbc_north_e 
    26632663 
    2664       SUBROUTINE mpp_lnk_bdy_3d( ptab, cd_type, psgn, ib_bdy ) 
     2664 
     2665   SUBROUTINE mpp_lnk_bdy_3d( ptab, cd_type, psgn, ib_bdy ) 
    26652666      !!---------------------------------------------------------------------- 
    26662667      !!                  ***  routine mpp_lnk_bdy_3d  *** 
     
    26832684      !! 
    26842685      !!---------------------------------------------------------------------- 
    2685  
    2686       USE lbcnfd          ! north fold 
    2687  
    2688       INCLUDE 'mpif.h' 
    2689  
    26902686      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied 
    26912687      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points 
     
    26942690      !                                                             ! =  1. , the sign is kept 
    26952691      INTEGER                         , INTENT(in   ) ::   ib_bdy   ! BDY boundary set 
     2692      ! 
    26962693      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices 
    2697       INTEGER  ::   imigr, iihom, ijhom        ! temporary integers 
     2694      INTEGER  ::   imigr, iihom, ijhom        ! local integers 
    26982695      INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend 
    2699       REAL(wp) ::   zland 
     2696      REAL(wp) ::   zland                      ! local scalar 
    27002697      INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend 
    27012698      ! 
    27022699      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ns, zt3sn   ! 3d for north-south & south-north 
    27032700      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ew, zt3we   ! 3d for east-west & west-east 
    2704  
    2705       !!---------------------------------------------------------------------- 
    2706        
     2701      !!---------------------------------------------------------------------- 
     2702      ! 
    27072703      ALLOCATE( zt3ns(jpi,jprecj,jpk,2), zt3sn(jpi,jprecj,jpk,2),   & 
    27082704         &      zt3ew(jpj,jpreci,jpk,2), zt3we(jpj,jpreci,jpk,2)  ) 
    27092705 
    2710       zland = 0.e0 
     2706      zland = 0.-WP 
    27112707 
    27122708      ! 1. standard boundary treatment 
    27132709      ! ------------------------------ 
    2714        
    27152710      !                                   ! East-West boundaries 
    27162711      !                                        !* Cyclic east-west 
    2717  
    27182712      IF( nbondi == 2) THEN 
    2719         IF (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN 
    2720           ptab( 1 ,:,:) = ptab(jpim1,:,:) 
    2721           ptab(jpi,:,:) = ptab(  2  ,:,:) 
    2722         ELSE 
    2723           IF( .NOT. cd_type == 'F' )   ptab(     1       :jpreci,:,:) = zland    ! south except F-point 
    2724           ptab(nlci-jpreci+1:jpi   ,:,:) = zland    ! north 
    2725         ENDIF 
     2713         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 
     2714            ptab( 1 ,:,:) = ptab(jpim1,:,:) 
     2715            ptab(jpi,:,:) = ptab(  2  ,:,:) 
     2716         ELSE 
     2717            IF( .NOT. cd_type == 'F' )   ptab(1:jpreci,:,:) = zland    ! south except F-point 
     2718            ptab(nlci-jpreci+1:jpi,:,:) = zland    ! north 
     2719         ENDIF 
    27262720      ELSEIF(nbondi == -1) THEN 
    2727         IF( .NOT. cd_type == 'F' )   ptab(     1       :jpreci,:,:) = zland    ! south except F-point 
     2721         IF( .NOT. cd_type == 'F' )   ptab(1:jpreci,:,:) = zland    ! south except F-point 
    27282722      ELSEIF(nbondi == 1) THEN 
    2729         ptab(nlci-jpreci+1:jpi   ,:,:) = zland    ! north 
     2723         ptab(nlci-jpreci+1:jpi,:,:) = zland    ! north 
    27302724      ENDIF                                     !* closed 
    27312725 
    27322726      IF (nbondj == 2 .OR. nbondj == -1) THEN 
    2733         IF( .NOT. cd_type == 'F' )   ptab(:,     1       :jprecj,:) = zland       ! south except F-point 
     2727        IF( .NOT. cd_type == 'F' )   ptab(:,1:jprecj,:) = zland       ! south except F-point 
    27342728      ELSEIF (nbondj == 2 .OR. nbondj == 1) THEN 
    2735         ptab(:,nlcj-jprecj+1:jpj   ,:) = zland       ! north 
    2736       ENDIF 
    2737        
    2738       ! 
    2739  
     2729        ptab(:,nlcj-jprecj+1:jpj,:) = zland       ! north 
     2730      ENDIF 
     2731      ! 
    27402732      ! 2. East and west directions exchange 
    27412733      ! ------------------------------------ 
     
    27942786      CASE ( 0 ) 
    27952787         DO jl = 1, jpreci 
    2796             ptab(jl      ,:,:) = zt3we(:,jl,:,2) 
     2788            ptab(      jl,:,:) = zt3we(:,jl,:,2) 
    27972789            ptab(iihom+jl,:,:) = zt3ew(:,jl,:,2) 
    27982790         END DO 
    27992791      CASE ( 1 ) 
    28002792         DO jl = 1, jpreci 
    2801             ptab(jl      ,:,:) = zt3we(:,jl,:,2) 
     2793            ptab(      jl,:,:) = zt3we(:,jl,:,2) 
    28022794         END DO 
    28032795      END SELECT 
     
    28852877   END SUBROUTINE mpp_lnk_bdy_3d 
    28862878 
    2887       SUBROUTINE mpp_lnk_bdy_2d( ptab, cd_type, psgn, ib_bdy ) 
     2879 
     2880   SUBROUTINE mpp_lnk_bdy_2d( ptab, cd_type, psgn, ib_bdy ) 
    28882881      !!---------------------------------------------------------------------- 
    28892882      !!                  ***  routine mpp_lnk_bdy_2d  *** 
     
    29062899      !! 
    29072900      !!---------------------------------------------------------------------- 
    2908  
    2909       USE lbcnfd          ! north fold 
    2910  
    2911       INCLUDE 'mpif.h' 
    2912  
    2913       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied 
    2914       CHARACTER(len=1)                , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points 
    2915       !                                                             ! = T , U , V , F , W points 
    2916       REAL(wp)                        , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary 
    2917       !                                                             ! =  1. , the sign is kept 
    2918       INTEGER                         , INTENT(in   ) ::   ib_bdy   ! BDY boundary set 
     2901      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied 
     2902      CHARACTER(len=1)            , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points 
     2903      !                                                         ! = T , U , V , F , W points 
     2904      REAL(wp)                    , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary 
     2905      !                                                         ! =  1. , the sign is kept 
     2906      INTEGER                     , INTENT(in   ) ::   ib_bdy   ! BDY boundary set 
     2907      ! 
    29192908      INTEGER  ::   ji, jj, jl             ! dummy loop indices 
    2920       INTEGER  ::   imigr, iihom, ijhom        ! temporary integers 
     2909      INTEGER  ::   imigr, iihom, ijhom        ! local integers 
    29212910      INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend 
    29222911      REAL(wp) ::   zland 
     
    29252914      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ns, zt2sn   ! 2d for north-south & south-north 
    29262915      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ew, zt2we   ! 2d for east-west & west-east 
    2927  
    29282916      !!---------------------------------------------------------------------- 
    29292917 
     
    29312919         &      zt2ew(jpj,jpreci,2), zt2we(jpj,jpreci,2)   ) 
    29322920 
    2933       zland = 0.e0 
     2921      zland = 0._wp 
    29342922 
    29352923      ! 1. standard boundary treatment 
    29362924      ! ------------------------------ 
    2937        
    29382925      !                                   ! East-West boundaries 
    2939       !                                        !* Cyclic east-west 
    2940  
    2941       IF( nbondi == 2) THEN 
    2942         IF (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN 
    2943           ptab( 1 ,:) = ptab(jpim1,:) 
    2944           ptab(jpi,:) = ptab(  2  ,:) 
    2945         ELSE 
    2946           IF( .NOT. cd_type == 'F' )   ptab(     1       :jpreci,:) = zland    ! south except F-point 
    2947           ptab(nlci-jpreci+1:jpi   ,:) = zland    ! north 
    2948         ENDIF 
     2926      !                                      !* Cyclic east-west 
     2927      IF( nbondi == 2 ) THEN 
     2928         IF (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN 
     2929            ptab( 1 ,:) = ptab(jpim1,:) 
     2930            ptab(jpi,:) = ptab(  2  ,:) 
     2931         ELSE 
     2932            IF(.NOT.cd_type == 'F' )  ptab(     1       :jpreci,:) = zland    ! south except F-point 
     2933                                      ptab(nlci-jpreci+1:jpi   ,:) = zland    ! north 
     2934         ENDIF 
    29492935      ELSEIF(nbondi == -1) THEN 
    2950         IF( .NOT. cd_type == 'F' )   ptab(     1       :jpreci,:) = zland    ! south except F-point 
     2936         IF( .NOT.cd_type == 'F' )    ptab(     1       :jpreci,:) = zland    ! south except F-point 
    29512937      ELSEIF(nbondi == 1) THEN 
    2952         ptab(nlci-jpreci+1:jpi   ,:) = zland    ! north 
    2953       ENDIF                                     !* closed 
    2954  
    2955       IF (nbondj == 2 .OR. nbondj == -1) THEN 
    2956         IF( .NOT. cd_type == 'F' )   ptab(:,     1       :jprecj) = zland       ! south except F-point 
     2938                                      ptab(nlci-jpreci+1:jpi   ,:) = zland    ! north 
     2939      ENDIF 
     2940      !                                      !* closed 
     2941      IF( nbondj == 2 .OR. nbondj == -1 ) THEN 
     2942         IF( .NOT.cd_type == 'F' )    ptab(:,     1       :jprecj) = zland    ! south except F-point 
    29572943      ELSEIF (nbondj == 2 .OR. nbondj == 1) THEN 
    2958         ptab(:,nlcj-jprecj+1:jpj) = zland       ! north 
    2959       ENDIF 
    2960        
    2961       ! 
    2962  
     2944                                      ptab(:,nlcj-jprecj+1:jpj   ) = zland    ! north 
     2945      ENDIF 
     2946      ! 
    29632947      ! 2. East and west directions exchange 
    29642948      ! ------------------------------------ 
     
    31073091      ! 
    31083092   END SUBROUTINE mpp_lnk_bdy_2d 
     3093 
    31093094 
    31103095   SUBROUTINE mpi_init_opa( ldtxt, ksft, code ) 
     
    31963181   END SUBROUTINE DDPDD_MPI 
    31973182 
     3183 
    31983184   SUBROUTINE mpp_lbc_north_icb( pt2d, cd_type, psgn, pr2dj) 
    31993185      !!--------------------------------------------------------------------- 
     
    32183204      !!                                                    ! north fold, =  1. otherwise 
    32193205      INTEGER, OPTIONAL       , INTENT(in   ) ::   pr2dj 
     3206      ! 
    32203207      INTEGER ::   ji, jj, jr 
    32213208      INTEGER ::   ierr, itaille, ildi, ilei, iilb 
     
    32243211      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE  ::  ztab_e, znorthloc_e 
    32253212      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE  ::  znorthgloio_e 
    3226  
    32273213      !!---------------------------------------------------------------------- 
    32283214      ! 
     
    32343220      ENDIF 
    32353221      ALLOCATE( ztab_e(jpiglo,4+2*ipr2dj), znorthloc_e(jpi,4+2*ipr2dj), znorthgloio_e(jpi,4+2*ipr2dj,jpni) ) 
    3236  
    3237       ! 
    3238       ztab_e(:,:) = 0.e0 
    3239  
    3240       ij=0 
     3222      ! 
     3223      ztab_e(:,:) = 0._wp 
     3224      ! 
     3225      ij = 0 
    32413226      ! put in znorthloc_e the last 4 jlines of pt2d 
    32423227      DO jj = nlcj - ijpj + 1 - ipr2dj, nlcj +ipr2dj 
     
    32803265      ! 
    32813266   END SUBROUTINE mpp_lbc_north_icb 
     3267 
    32823268 
    32833269   SUBROUTINE mpp_lnk_2d_icb( pt2d, cd_type, psgn, jpri, jprj ) 
     
    33003286      !!                    noso   : number for local neighboring processors 
    33013287      !!                    nono   : number for local neighboring processors 
    3302       !! 
    33033288      !!---------------------------------------------------------------------- 
    33043289      INTEGER                                             , INTENT(in   ) ::   jpri 
     
    34593444 
    34603445   END SUBROUTINE mpp_lnk_2d_icb 
     3446    
    34613447#else 
    34623448   !!---------------------------------------------------------------------- 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

    r5836 r5883  
    44   !! Ocean forcing:  read input field for surface boundary condition 
    55   !!===================================================================== 
    6    !! History :  2.0  !  06-2006  (S. Masson, G. Madec) Original code 
    7    !!                 !  05-2008  (S. Alderson) Modified for Interpolation in memory 
    8    !!                 !                         from input grid to model grid 
    9    !!                 !  10-2013  (D. Delrosso, P. Oddo) implement suppression of  
    10    !!                 !                         land point prior to interpolation 
     6   !! History :  2.0  !  06-2006  (S. Masson, G. Madec)  Original code 
     7   !!                 !  05-2008  (S. Alderson)  Modified for Interpolation in memory from input grid to model grid 
     8   !!                 !  10-2013  (D. Delrosso, P. Oddo)  suppression of land point prior to interpolation 
    119   !!---------------------------------------------------------------------- 
    1210 
     
    1513   !!                   surface boundary condition 
    1614   !!---------------------------------------------------------------------- 
    17    USE oce             ! ocean dynamics and tracers 
    18    USE dom_oce         ! ocean space and time domain 
    19    USE phycst          ! ??? 
    20    USE in_out_manager  ! I/O manager 
    21    USE iom             ! I/O manager library 
    22    USE geo2ocean       ! for vector rotation on to model grid 
    23    USE lib_mpp         ! MPP library 
    24    USE wrk_nemo        ! work arrays 
    25    USE lbclnk          ! ocean lateral boundary conditions (C1D case) 
    26    USE ioipsl, ONLY ymds2ju, ju2ymds   ! for calendar 
     15   USE oce            ! ocean dynamics and tracers 
     16   USE dom_oce        ! ocean space and time domain 
     17   USE phycst         ! physical constant 
     18   USE in_out_manager ! I/O manager 
     19   USE iom            ! I/O manager library 
     20   USE geo2ocean      ! for vector rotation on to model grid 
     21   USE lib_mpp        ! MPP library 
     22   USE wrk_nemo       ! work arrays 
     23   USE lbclnk         ! ocean lateral boundary conditions (C1D case) 
     24   USE ioipsl, ONLY   : ymds2ju, ju2ymds   ! for calendar 
    2725   USE sbc_oce 
    2826    
     
    6058      INTEGER , DIMENSION(2)          ::   nrec_b       ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year) 
    6159      INTEGER , DIMENSION(2)          ::   nrec_a       ! after  record (1: index, 2: second since Jan. 1st 00h of nit000 year) 
    62       REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:  ) ::   fnow       ! input fields interpolated to now time step 
    63       REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta       ! 2 consecutive record of input fields 
     60      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:  ) ::   fnow   ! input fields interpolated to now time step 
     61      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta   ! 2 consecutive record of input fields 
    6462      CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key 
    6563      !                                                 ! into the WGTLIST structure 
     
    133131      INTEGER  , INTENT(in   ), OPTIONAL     ::   kit       ! subcycle timestep for timesplitting option 
    134132      INTEGER  , INTENT(in   ), OPTIONAL     ::   kt_offset ! provide fields at time other than "now" 
    135                                                             !   kt_offset = -1 => fields at "before" time level 
    136                                                             !   kt_offset = +1 => fields at "after"  time level 
    137                                                             !   etc. 
    138       !! 
    139       INTEGER  ::   itmp       ! temporary variable 
    140       INTEGER  ::   imf        ! size of the structure sd 
    141       INTEGER  ::   jf         ! dummy indices 
    142       INTEGER  ::   isecend    ! number of second since Jan. 1st 00h of nit000 year at nitend 
    143       INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step 
    144       INTEGER  ::   it_offset  ! local time offset variable 
    145       LOGICAL  ::   llnxtyr    ! open next year  file? 
    146       LOGICAL  ::   llnxtmth   ! open next month file? 
    147       LOGICAL  ::   llstop     ! stop is the file does not exist 
     133      !                                                     !   kt_offset = -1 => fields at "before" time level 
     134      !                                                     !   kt_offset = +1 => fields at "after"  time level 
     135      !                                                     !   etc. 
     136      ! 
     137      INTEGER  ::   itmp         ! temporary variable 
     138      INTEGER  ::   imf          ! size of the structure sd 
     139      INTEGER  ::   jf           ! dummy indices 
     140      INTEGER  ::   isecend      ! number of second since Jan. 1st 00h of nit000 year at nitend 
     141      INTEGER  ::   isecsbc      ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step 
     142      INTEGER  ::   it_offset    ! local time offset variable 
     143      LOGICAL  ::   llnxtyr      ! open next year  file? 
     144      LOGICAL  ::   llnxtmth     ! open next month file? 
     145      LOGICAL  ::   llstop       ! stop is the file does not exist 
    148146      LOGICAL  ::   ll_firstcall ! true if this is the first call to fld_read for this set of fields 
    149       REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation 
    150       REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation 
    151       CHARACTER(LEN=1000) ::   clfmt   ! write format 
    152       TYPE(MAP_POINTER) ::   imap   ! global-to-local mapping indices 
     147      REAL(wp) ::   ztinta       ! ratio applied to after  records when doing time interpolation 
     148      REAL(wp) ::   ztintb       ! ratio applied to before records when doing time interpolation 
     149      CHARACTER(LEN=1000) ::   clfmt  ! write format 
     150      TYPE(MAP_POINTER)   ::   imap   ! global-to-local mapping indices 
    153151      !!--------------------------------------------------------------------- 
    154152      ll_firstcall = kt == nit000 
     
    299297         END DO                                    ! --- end loop over field --- ! 
    300298         ! 
    301          !                                         ! ====================================== ! 
    302       ENDIF                                        ! update field at each kn_fsbc time-step ! 
    303       !                                            ! ====================================== ! 
     299      ENDIF 
    304300      ! 
    305301   END SUBROUTINE fld_read 
     
    333329      llprevday  = .FALSE. 
    334330      isec_week  = 0 
    335              
     331      ! 
    336332      ! define record informations 
    337333      CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. )  ! return before values in sdjf%nrec_a (as we will swap it later) 
    338  
     334      ! 
    339335      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar  
    340  
     336      ! 
    341337      IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure 
    342  
     338         ! 
    343339         IF( sdjf%nrec_a(1) == 0  ) THEN   ! we redefine record sdjf%nrec_a(1) with the last record of previous year file 
    344340            IF    ( sdjf%nfreqh == -12 ) THEN   ! yearly mean 
     
    391387         ! 
    392388         CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev ) 
    393  
     389         ! 
    394390         ! if previous year/month/day file does not exist, we switch to the current year/month/day 
    395391         IF( llprev .AND. sdjf%num <= 0 ) THEN 
     
    401397            CALL fld_clopn( sdjf ) 
    402398         ENDIF 
    403           
     399         ! 
    404400         IF( llprev ) THEN   ! check if the record sdjf%nrec_a(1) exists in the file 
    405401            idvar = iom_varid( sdjf%num, sdjf%clvar )                                        ! id of the variable sdjf%clvar 
     
    408404            sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec )   ! make sure we select an existing record 
    409405         ENDIF 
    410  
    411          ! read before data in after arrays(as we will swap it later) 
    412          CALL fld_get( sdjf, map ) 
    413  
     406         ! 
     407         CALL fld_get( sdjf, map )         ! read before data in after arrays(as we will swap it later) 
     408         ! 
    414409         clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" 
    415410         IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday 
    416  
     411         ! 
    417412      ENDIF 
    418413      ! 
     
    435430      LOGICAL  , INTENT(in   ), OPTIONAL ::   ldbefore  ! sent back before record values (default = .FALSE.) 
    436431      INTEGER  , INTENT(in   ), OPTIONAL ::   kit       ! index of barotropic subcycle 
    437                                                         ! used only if sdjf%ln_tint = .TRUE. 
     432      !                                                 ! used only if sdjf%ln_tint = .TRUE. 
    438433      INTEGER  , INTENT(in   ), OPTIONAL ::   kt_offset ! Offset of required time level compared to "now" 
    439                                                         !   time level in units of time steps. 
    440       !! 
     434      !                                                 !   time level in units of time steps. 
     435      ! 
    441436      LOGICAL  ::   llbefore    ! local definition of ldbefore 
    442437      INTEGER  ::   iendrec     ! end of this record (in seconds) 
     
    592587      !! ** Purpose :   read the data 
    593588      !!---------------------------------------------------------------------- 
    594       TYPE(FLD), INTENT(inout) ::   sdjf   ! input field related variables 
    595       TYPE(MAP_POINTER),INTENT(in) ::   map   ! global-to-local mapping indices 
    596       !! 
    597       INTEGER                  ::   ipk    ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
    598       INTEGER                  ::   iw     ! index into wgts array 
    599       INTEGER                  ::   ipdom  ! index of the domain 
    600       INTEGER                  ::   idvar  ! variable ID 
    601       INTEGER                  ::   idmspc ! number of spatial dimensions 
    602       LOGICAL                  ::   lmoor  ! C1D case: point data 
     589      TYPE(FLD)        , INTENT(inout) ::   sdjf   ! input field related variables 
     590      TYPE(MAP_POINTER), INTENT(in   ) ::   map    ! global-to-local mapping indices 
     591      ! 
     592      INTEGER ::   ipk      ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
     593      INTEGER ::   iw       ! index into wgts array 
     594      INTEGER ::   ipdom    ! index of the domain 
     595      INTEGER ::   idvar    ! variable ID 
     596      INTEGER ::   idmspc  ! number of spatial dimensions 
     597      LOGICAL ::   lmoor    ! C1D case: point data 
    603598      !!--------------------------------------------------------------------- 
    604599      ! 
     
    611606      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 
    612607         CALL wgt_list( sdjf, iw ) 
    613          IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fdta(:,:,:,2), &  
    614               & sdjf%nrec_a(1), sdjf%lsmname ) 
    615          ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fnow(:,:,:  ), & 
    616               & sdjf%nrec_a(1), sdjf%lsmname ) 
     608         IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fdta(:,:,:,2),          &  
     609            &                                                                          sdjf%nrec_a(1), sdjf%lsmname ) 
     610         ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fnow(:,:,:  ),          & 
     611            &                                                                          sdjf%nrec_a(1), sdjf%lsmname ) 
    617612         ENDIF 
    618613      ELSE 
    619          IF( SIZE(sdjf%fnow, 1) == jpi ) THEN  ;  ipdom = jpdom_data 
    620          ELSE                                  ;  ipdom = jpdom_unknown 
     614         IF( SIZE(sdjf%fnow, 1) == jpi ) THEN   ;   ipdom = jpdom_data 
     615         ELSE                                   ;   ipdom = jpdom_unknown 
    621616         ENDIF 
    622617         ! C1D case: If product of spatial dimensions == ipk, then x,y are of 
    623618         ! size 1 (point/mooring data): this must be read onto the central grid point 
    624619         idvar  = iom_varid( sdjf%num, sdjf%clvar ) 
    625          idmspc = iom_file( sdjf%num )%ndims( idvar ) 
     620         idmspc = iom_file ( sdjf%num )%ndims( idvar ) 
    626621         IF( iom_file( sdjf%num )%luld( idvar ) )   idmspc = idmspc - 1 
    627          lmoor  = (idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk) 
     622         lmoor  = (  idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk  ) 
    628623         ! 
    629624         SELECT CASE( ipk ) 
     
    660655      ! 
    661656      sdjf%rotn(2) = .false.   ! vector not yet rotated 
    662  
     657      ! 
    663658   END SUBROUTINE fld_get 
     659 
    664660 
    665661   SUBROUTINE fld_map( num, clvar, dta, nrec, map ) 
     
    688684      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read  ! work space for global data 
    689685      !!--------------------------------------------------------------------- 
    690              
     686      ! 
    691687      ipi = SIZE( dta, 1 ) 
    692688      ipj = 1 
    693689      ipk = SIZE( dta, 3 ) 
    694  
     690      ! 
    695691      idvar   = iom_varid( num, clvar ) 
    696692      ilendta = iom_file(num)%dimsz(1,idvar) 
     
    698694#if defined key_bdy 
    699695      ipj = iom_file(num)%dimsz(2,idvar) 
    700       IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 
     696      IF( map%ll_unstruc) THEN  ! unstructured open boundary data file 
    701697         dta_read => dta_global 
    702       ELSE                      ! structured open boundary data file 
     698      ELSE                       ! structured open boundary data file 
    703699         dta_read => dta_global2 
    704700      ENDIF 
    705701#endif 
    706702 
    707       IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta 
     703      IF(lwp) WRITE(numout,*) 'Dim size for ',        TRIM(clvar),' is ', ilendta 
    708704      IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk 
    709  
     705      ! 
    710706      SELECT CASE( ipk ) 
    711707      CASE(1)        ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec ) 
     
    713709      END SELECT 
    714710      ! 
    715       IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file 
     711      IF( map%ll_unstruc ) THEN ! unstructured open boundary data file 
    716712         DO ib = 1, ipi 
    717713            DO ik = 1, ipk 
     
    728724         END DO 
    729725      ENDIF 
    730  
     726      ! 
    731727   END SUBROUTINE fld_map 
    732728 
     
    738734      !! ** Purpose :   Vector fields may need to be rotated onto the local grid direction 
    739735      !!---------------------------------------------------------------------- 
    740       INTEGER  , INTENT(in   )               ::   kt        ! ocean time step 
    741       TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables 
    742       !! 
    743       INTEGER                           ::   ju,jv,jk,jn  ! loop indices 
    744       INTEGER                           ::   imf          ! size of the structure sd 
    745       INTEGER                           ::   ill          ! character length 
    746       INTEGER                           ::   iv           ! indice of V component 
     736      INTEGER                , INTENT(in   ) ::   kt   ! ocean time step 
     737      TYPE(FLD), DIMENSION(:), INTENT(inout) ::   sd   ! input field related variables 
     738      ! 
     739      INTEGER ::   ju, jv, jk, jn  ! loop indices 
     740      INTEGER ::   imf             ! size of the structure sd 
     741      INTEGER ::   ill             ! character length 
     742      INTEGER ::   iv              ! indice of V component 
     743      CHARACTER (LEN=100)               ::   clcomp       ! dummy weight name 
    747744      REAL(wp), POINTER, DIMENSION(:,:) ::   utmp, vtmp   ! temporary arrays for vector rotation 
    748       CHARACTER (LEN=100)               ::   clcomp       ! dummy weight name 
    749       !!--------------------------------------------------------------------- 
    750  
    751       CALL wrk_alloc( jpi,jpj, utmp, vtmp ) 
    752  
     745      !!--------------------------------------------------------------------- 
     746      ! 
     747      CALL wrk_alloc( jpi,jpj,   utmp, vtmp ) 
     748      ! 
    753749      !! (sga: following code should be modified so that pairs arent searched for each time 
    754750      ! 
     
    786782       END DO 
    787783      ! 
    788       CALL wrk_dealloc( jpi,jpj, utmp, vtmp ) 
     784      CALL wrk_dealloc( jpi,jpj,   utmp, vtmp ) 
    789785      ! 
    790786   END SUBROUTINE fld_rot 
     
    802798      INTEGER, OPTIONAL, INTENT(in   ) ::   kday     ! day value 
    803799      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.) 
    804       !! 
     800      ! 
    805801      LOGICAL :: llprevyr              ! are we reading previous year  file? 
    806802      LOGICAL :: llprevmth             ! are we reading previous month file? 
     
    853849      ! 
    854850      IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN   ! new file to be open  
    855  
     851         ! 
    856852         sdjf%clname = TRIM(clname) 
    857853         IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open 
    858854         CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 ) 
    859  
     855         ! 
    860856         ! find the last record to be read -> update sdjf%nreclast 
    861857         indexyr = iyear - nyear + 1 
     
    866862         CASE ( 2 )   ;   imonth_len = 31   ! next     year -> imonth = 1 
    867863         END SELECT 
    868           
     864         ! 
    869865         ! last record to be read in the current file 
    870866         IF    ( sdjf%nfreqh == -12 ) THEN                 ;   sdjf%nreclast = 1    !  yearly mean 
     
    880876            ENDIF 
    881877         ENDIF 
    882           
     878         ! 
    883879      ENDIF 
    884880      ! 
     
    901897      INTEGER  ::   jf       ! dummy indices 
    902898      !!--------------------------------------------------------------------- 
    903  
     899      ! 
    904900      DO jf = 1, SIZE(sdf) 
    905901         sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname ) 
     
    923919         sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 
    924920      END DO 
    925  
     921      ! 
    926922      IF(lwp) THEN      ! control print 
    927923         WRITE(numout,*) 
     
    943939         END DO 
    944940      ENDIF 
    945        
     941      ! 
    946942   END SUBROUTINE fld_fill 
    947943 
     
    958954      TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file 
    959955      INTEGER    , INTENT(inout) ::   kwgt      ! index of weights 
    960       !! 
     956      ! 
    961957      INTEGER ::   kw, nestid   ! local integer 
    962958      LOGICAL ::   found        ! local logical 
     
    966962      !! weights filename is either present or we hit the end of the list 
    967963      found = .FALSE. 
    968  
     964      ! 
    969965      !! because agrif nest part of filenames are now added in iom_open 
    970966      !! to distinguish between weights files on the different grids, need to track 
     
    10281024      TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file 
    10291025      !! 
    1030       INTEGER                           ::   jn            ! dummy loop indices 
    1031       INTEGER                           ::   inum          ! temporary logical unit 
    1032       INTEGER                           ::   id            ! temporary variable id 
    1033       INTEGER                           ::   ipk           ! temporary vertical dimension 
    1034       CHARACTER (len=5)                 ::   aname 
     1026      INTEGER ::   jn         ! dummy loop indices 
     1027      INTEGER ::   inum       ! local logical unit 
     1028      INTEGER ::   id         ! local variable id 
     1029      INTEGER ::   ipk        ! local vertical dimension 
     1030      INTEGER ::   zwrap      ! local integer 
     1031      LOGICAL ::   cyclical   !  
     1032      CHARACTER (len=5) ::   aname   ! 
    10351033      INTEGER , DIMENSION(:), ALLOCATABLE ::   ddims 
    10361034      INTEGER , POINTER, DIMENSION(:,:) ::   data_src 
    10371035      REAL(wp), POINTER, DIMENSION(:,:) ::   data_tmp 
    1038       LOGICAL                           ::   cyclical 
    1039       INTEGER                           ::   zwrap      ! local integer 
    1040       !!---------------------------------------------------------------------- 
    1041       ! 
    1042       CALL wrk_alloc( jpi,jpj, data_src )   ! integer 
    1043       CALL wrk_alloc( jpi,jpj, data_tmp ) 
     1036      !!---------------------------------------------------------------------- 
     1037      ! 
     1038      CALL wrk_alloc( jpi,jpj,   data_src )   ! integer 
     1039      CALL wrk_alloc( jpi,jpj,   data_tmp ) 
    10441040      ! 
    10451041      IF( nxt_wgt > tot_wgts ) THEN 
     
    11511147         ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) ) 
    11521148         IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) ) 
    1153  
     1149         ! 
    11541150         nxt_wgt = nxt_wgt + 1 
    1155  
     1151         ! 
    11561152      ELSE  
    11571153         CALL ctl_stop( '    fld_weight : unable to read the file ' ) 
     
    11661162 
    11671163 
    1168    SUBROUTINE apply_seaoverland(clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm, & 
    1169                           &      jpj2_lsm,itmpi,itmpj,itmpz,rec1_lsm,recn_lsm) 
     1164   SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm,  & 
     1165                          &      jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm ) 
    11701166      !!--------------------------------------------------------------------- 
    11711167      !!                    ***  ROUTINE apply_seaoverland  *** 
     
    11761172      !!      D. Delrosso INGV           
    11771173      !!----------------------------------------------------------------------  
    1178       INTEGER                                   :: inum,jni,jnj,jnz,jc                  ! temporary indices 
    1179       INTEGER,                   INTENT(in)     :: itmpi,itmpj,itmpz                    ! lengths 
    1180       INTEGER,                   INTENT(in)     :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm  ! temporary indices 
    1181       INTEGER, DIMENSION(3),     INTENT(in)     :: rec1_lsm,recn_lsm                    ! temporary arrays for start and length 
    1182       REAL(wp),DIMENSION (:,:,:),INTENT(inout)  :: zfieldo                              ! input/output array for seaoverland application 
    1183       REAL(wp),DIMENSION (:,:,:),ALLOCATABLE    :: zslmec1                              ! temporary array for land point detection 
    1184       REAL(wp),DIMENSION (:,:),  ALLOCATABLE    :: zfieldn                              ! array of forcing field with undeff for land points 
    1185       REAL(wp),DIMENSION (:,:),  ALLOCATABLE    :: zfield                               ! array of forcing field 
    1186       CHARACTER (len=100),       INTENT(in)     :: clmaskfile                           ! land/sea mask file name 
    1187       !!--------------------------------------------------------------------- 
    1188       ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) ) 
    1189       ALLOCATE ( zfieldn(itmpi,itmpj) ) 
    1190       ALLOCATE ( zfield(itmpi,itmpj) ) 
    1191  
     1174      INTEGER,                   INTENT(in   ) :: itmpi,itmpj,itmpz                    ! lengths 
     1175      INTEGER,                   INTENT(in   ) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm  ! temporary indices 
     1176      INTEGER, DIMENSION(3),     INTENT(in   ) :: rec1_lsm,recn_lsm                    ! temporary arrays for start and length 
     1177      REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo                              ! input/output array for seaoverland application 
     1178      CHARACTER (len=100),       INTENT(in   ) :: clmaskfile                           ! land/sea mask file name 
     1179      ! 
     1180      INTEGER :: inum,jni,jnj,jnz,jc   ! local indices 
     1181      REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1             ! local array for land point detection 
     1182      REAL(wp),DIMENSION (:,:),  ALLOCATABLE :: zfieldn   ! array of forcing field with undeff for land points 
     1183      REAL(wp),DIMENSION (:,:),  ALLOCATABLE :: zfield    ! array of forcing field 
     1184      !!--------------------------------------------------------------------- 
     1185      ! 
     1186      ALLOCATE ( zslmec1(itmpi,itmpj,itmpz), zfieldn(itmpi,itmpj), zfield(itmpi,itmpj) ) 
     1187      ! 
    11921188      ! Retrieve the land sea mask data 
    11931189      CALL iom_open( clmaskfile, inum ) 
    11941190      SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) ) 
    11951191      CASE(1) 
    1196            CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm) 
     1192         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm) 
    11971193      CASE DEFAULT 
    1198            CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm) 
     1194         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm) 
    11991195      END SELECT 
    12001196      CALL iom_close( inum ) 
    1201  
    1202       DO jnz=1,rec1_lsm(3)                            !! Loop over k dimension 
    1203  
    1204          DO jni=1,itmpi                               !! copy the original field into a tmp array 
    1205             DO jnj=1,itmpj                            !! substituting undeff over land points 
    1206             zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz) 
    1207                IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN 
    1208                   zfieldn(jni,jnj) = undeff_lsm 
    1209                ENDIF 
     1197      ! 
     1198      DO jnz=1,rec1_lsm(3)             !! Loop over k dimension 
     1199         ! 
     1200         DO jni = 1, itmpi                               !! copy the original field into a tmp array 
     1201            DO jnj = 1, itmpj                            !! substituting undeff over land points 
     1202               zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz) 
     1203               IF( zslmec1(jni,jnj,jnz) == 1. )   zfieldn(jni,jnj) = undeff_lsm 
    12101204            END DO 
    12111205         END DO 
    1212    
    1213       CALL seaoverland(zfieldn,itmpi,itmpj,zfield) 
    1214       DO jc=1,nn_lsm 
    1215          CALL seaoverland(zfield,itmpi,itmpj,zfield) 
    1216       END DO 
    1217  
    1218       !   Check for Undeff and substitute original values 
    1219       IF(ANY(zfield==undeff_lsm)) THEN 
    1220          DO jni=1,itmpi 
    1221             DO jnj=1,itmpj 
    1222                IF (zfield(jni,jnj)==undeff_lsm) THEN 
    1223                   zfield(jni,jnj) = zfieldo(jni,jnj,jnz) 
    1224                ENDIF 
    1225             ENDDO 
    1226          ENDDO 
    1227       ENDIF 
    1228  
    1229       zfieldo(:,:,jnz)=zfield(:,:) 
    1230  
    1231       END DO                          !! End Loop over k dimension 
    1232  
    1233       DEALLOCATE ( zslmec1 ) 
    1234       DEALLOCATE ( zfieldn ) 
    1235       DEALLOCATE ( zfield ) 
    1236  
     1206         ! 
     1207         CALL seaoverland( zfieldn, itmpi, itmpj, zfield ) 
     1208         DO jc = 1, nn_lsm 
     1209            CALL seaoverland( zfield, itmpi, itmpj, zfield ) 
     1210         END DO 
     1211         ! 
     1212         !   Check for Undeff and substitute original values 
     1213         IF( ANY(zfield==undeff_lsm) ) THEN 
     1214            DO jni = 1, itmpi 
     1215               DO jnj = 1, itmpj 
     1216                  IF( zfield(jni,jnj)==undeff_lsm )   zfield(jni,jnj) = zfieldo(jni,jnj,jnz) 
     1217               END DO 
     1218            END DO 
     1219         ENDIF 
     1220         ! 
     1221         zfieldo(:,:,jnz) = zfield(:,:) 
     1222         ! 
     1223      END DO                           !! End Loop over k dimension 
     1224      ! 
     1225      DEALLOCATE ( zslmec1, zfieldn, zfield ) 
     1226      ! 
    12371227   END SUBROUTINE apply_seaoverland  
    12381228 
    12391229 
    1240    SUBROUTINE seaoverland(zfieldn,ileni,ilenj,zfield) 
     1230   SUBROUTINE seaoverland( zfieldn, ileni, ilenj, zfield ) 
    12411231      !!--------------------------------------------------------------------- 
    12421232      !!                    ***  ROUTINE seaoverland  *** 
     
    12451235      !!      D. Delrosso INGV 
    12461236      !!----------------------------------------------------------------------  
    1247       INTEGER,INTENT(in)                       :: ileni,ilenj              ! lengths  
    1248       REAL,DIMENSION (ileni,ilenj),INTENT(in)  :: zfieldn                  ! array of forcing field with undeff for land points 
    1249       REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield                   ! array of forcing field 
    1250       REAL,DIMENSION (ileni,ilenj)             :: zmat1,zmat2,zmat3,zmat4  ! temporary arrays for seaoverland application 
    1251       REAL,DIMENSION (ileni,ilenj)             :: zmat5,zmat6,zmat7,zmat8  ! temporary arrays for seaoverland application 
    1252       REAL,DIMENSION (ileni,ilenj)             :: zlsm2d                   ! temporary arrays for seaoverland application 
    1253       REAL,DIMENSION (ileni,ilenj,8)           :: zlsm3d                   ! temporary arrays for seaoverland application 
    1254       LOGICAL,DIMENSION (ileni,ilenj,8)        :: ll_msknan3d              ! logical mask for undeff detection 
    1255       LOGICAL,DIMENSION (ileni,ilenj)          :: ll_msknan2d              ! logical mask for undeff detection 
     1237      INTEGER                      , INTENT(in   ) :: ileni,ilenj   ! lengths  
     1238      REAL, DIMENSION (ileni,ilenj), INTENT(in   ) :: zfieldn       ! array of forcing field with undeff for land points 
     1239      REAL, DIMENSION (ileni,ilenj), INTENT(  out) :: zfield        ! array of forcing field 
     1240      ! 
     1241      REAL   , DIMENSION (ileni,ilenj)   :: zmat1, zmat2, zmat3, zmat4  ! local arrays  
     1242      REAL   , DIMENSION (ileni,ilenj)   :: zmat5, zmat6, zmat7, zmat8  !   -     -  
     1243      REAL   , DIMENSION (ileni,ilenj)   :: zlsm2d                      !   -     -  
     1244      REAL   , DIMENSION (ileni,ilenj,8) :: zlsm3d                      !   -     - 
     1245      LOGICAL, DIMENSION (ileni,ilenj,8) :: ll_msknan3d                 ! logical mask for undeff detection 
     1246      LOGICAL, DIMENSION (ileni,ilenj)   :: ll_msknan2d                 ! logical mask for undeff detection 
    12561247      !!----------------------------------------------------------------------  
    1257       zmat8 = eoshift(zfieldn   ,  SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/)    ,DIM=2) 
    1258       zmat1 = eoshift(zmat8     ,  SHIFT=-1, BOUNDARY = (/zmat8(1,:)/)      ,DIM=1) 
    1259       zmat2 = eoshift(zfieldn   ,  SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/)    ,DIM=1) 
    1260       zmat4 = eoshift(zfieldn   ,  SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2) 
    1261       zmat3 = eoshift(zmat4     ,  SHIFT=-1, BOUNDARY = (/zmat4(1,:)/)      ,DIM=1) 
    1262       zmat5 = eoshift(zmat4     ,  SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/)  ,DIM=1) 
    1263       zmat6 = eoshift(zfieldn   ,  SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1) 
    1264       zmat7 = eoshift(zmat8     ,  SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/)  ,DIM=1) 
    1265  
     1248      zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/)     , DIM=2 ) 
     1249      zmat1 = eoshift( zmat8   , SHIFT=-1 , BOUNDARY = (/zmat8(1,:)/)       , DIM=1 ) 
     1250      zmat2 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(1,:)/)     , DIM=1 ) 
     1251      zmat4 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(:,ilenj)/) , DIM=2 ) 
     1252      zmat3 = eoshift( zmat4   , SHIFT=-1 , BOUNDARY = (/zmat4(1,:)/)       , DIM=1 ) 
     1253      zmat5 = eoshift( zmat4   , SHIFT= 1 , BOUNDARY = (/zmat4(ileni,:)/)   , DIM=1 ) 
     1254      zmat6 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(ileni,:)/) , DIM=1 ) 
     1255      zmat7 = eoshift( zmat8   , SHIFT= 1 , BOUNDARY = (/zmat8(ileni,:)/)   , DIM=1 ) 
     1256      ! 
    12661257      zlsm3d  = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /)) 
    1267       ll_msknan3d = .not.(zlsm3d==undeff_lsm) 
    1268       ll_msknan2d = .not.(zfieldn==undeff_lsm)  ! FALSE where is Undeff (land) 
    1269       zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 ))   )) 
    1270       WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp)  zlsm2d = undeff_lsm 
    1271       zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d) 
     1258      ll_msknan3d = .NOT.( zlsm3d  == undeff_lsm ) 
     1259      ll_msknan2d = .NOT.( zfieldn == undeff_lsm )  ! FALSE where is Undeff (land) 
     1260      zlsm2d = SUM( zlsm3d, 3 , ll_msknan3d ) / MAX( 1 , COUNT( ll_msknan3d , 3 ) ) 
     1261      WHERE( COUNT( ll_msknan3d , 3 ) == 0._wp )   zlsm2d = undeff_lsm 
     1262      zfield = MERGE( zfieldn, zlsm2d, ll_msknan2d ) 
     1263      ! 
    12721264   END SUBROUTINE seaoverland 
    12731265 
     
    12881280      INTEGER                   , INTENT(in   ) ::   nrec    ! record number to read (ie time slice) 
    12891281      CHARACTER(LEN=*)          , INTENT(in   ) ::   lsmfile ! land sea mask file name 
    1290       !!  
    1291       REAL(wp),DIMENSION(:,:,:),ALLOCATABLE     ::   ztmp_fly_dta,zfieldo                  ! temporary array of values on input grid 
    1292       INTEGER, DIMENSION(3)                     ::   rec1,recn                             ! temporary arrays for start and length 
    1293       INTEGER, DIMENSION(3)                     ::   rec1_lsm,recn_lsm                     ! temporary arrays for start and length in case of seaoverland 
    1294       INTEGER                                   ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2       ! temporary indices 
    1295       INTEGER                                   ::   jk, jn, jm, jir, jjr                  ! loop counters 
    1296       INTEGER                                   ::   ni, nj                                ! lengths 
    1297       INTEGER                                   ::   jpimin,jpiwid                         ! temporary indices 
    1298       INTEGER                                   ::   jpimin_lsm,jpiwid_lsm                 ! temporary indices 
    1299       INTEGER                                   ::   jpjmin,jpjwid                         ! temporary indices 
    1300       INTEGER                                   ::   jpjmin_lsm,jpjwid_lsm                 ! temporary indices 
    1301       INTEGER                                   ::   jpi1,jpi2,jpj1,jpj2                   ! temporary indices 
    1302       INTEGER                                   ::   jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm   ! temporary indices 
    1303       INTEGER                                   ::   itmpi,itmpj,itmpz                     ! lengths 
    1304        
     1282      ! 
     1283      INTEGER, DIMENSION(3) ::   rec1, recn           ! temporary arrays for start and length 
     1284      INTEGER, DIMENSION(3) ::   rec1_lsm, recn_lsm   ! temporary arrays for start and length in case of seaoverland 
     1285      INTEGER ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2    ! temporary indices 
     1286      INTEGER ::   jk, jn, jm, jir, jjr               ! loop counters 
     1287      INTEGER ::   ni, nj                             ! lengths 
     1288      INTEGER ::   jpimin,jpiwid                      ! temporary indices 
     1289      INTEGER ::   jpimin_lsm,jpiwid_lsm              ! temporary indices 
     1290      INTEGER ::   jpjmin,jpjwid                      ! temporary indices 
     1291      INTEGER ::   jpjmin_lsm,jpjwid_lsm              ! temporary indices 
     1292      INTEGER ::   jpi1,jpi2,jpj1,jpj2                ! temporary indices 
     1293      INTEGER ::   jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm   ! temporary indices 
     1294      INTEGER ::   itmpi,itmpj,itmpz                     ! lengths 
     1295      REAL(wp),DIMENSION(:,:,:), ALLOCATABLE ::   ztmp_fly_dta, zfieldo                 ! local array of values on input grid      
    13051296      !!---------------------------------------------------------------------- 
    13061297      ! 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90

    r5836 r5883  
    55   !!====================================================================== 
    66   !! History :  OPA  !  07-1996  (O. Marti)  Original code 
    7    !!   NEMO     1.0  !  02-2008  (G. Madec)  F90: Free form 
    8    !!            3.0  !   
     7   !!   NEMO     1.0  !  06-2006  (G. Madec )  Free form, F90 + opt. 
     8   !!                 !  04-2007  (S. Masson)  angle: Add T, F points and bugfix in cos lateral boundary 
     9   !!            3.0  !  07-2008  (G. Madec)  geo2oce suppress lon/lat agruments 
     10   !!            3.7  !  11-2015  (G. Madec)  remove the unused repere and repcmo routines 
    911   !!---------------------------------------------------------------------- 
    1012 
    1113   !!---------------------------------------------------------------------- 
    12    !!   repcmo      :  
    13    !!   angle       : 
    14    !!   geo2oce     : 
    15    !!   repere      :   old routine suppress it ??? 
     14   !!   rot_rep       : Rotate the Repere: geographic grid <==> stretched coordinates grid 
     15   !!   angle         : 
     16   !!   geo2oce       : 
     17   !!   oce2geo       : 
    1618   !!---------------------------------------------------------------------- 
    17    USE dom_oce         ! mesh and scale factors 
    18    USE phycst          ! physical constants 
    19    USE in_out_manager  ! I/O manager 
    20    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    21    USE lib_mpp         ! MPP library 
     19   USE dom_oce        ! mesh and scale factors 
     20   USE phycst         ! physical constants 
     21   ! 
     22   USE in_out_manager ! I/O manager 
     23   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     24   USE lib_mpp        ! MPP library 
    2225 
    2326   IMPLICIT NONE 
    2427   PRIVATE 
    2528 
    26    PUBLIC   rot_rep, repcmo, repere, geo2oce, oce2geo   ! only rot_rep should be used 
    27                                              ! repcmo and repere are keep only for compatibility. 
    28                                              ! they are only a useless overlay of rot_rep 
    29  
    30    PUBLIC   obs_rot 
    31  
    32    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   & 
    33       gsint, gcost,   &  ! cos/sin between model grid lines and NP direction at T point 
    34       gsinu, gcosu,   &  ! cos/sin between model grid lines and NP direction at U point 
    35       gsinv, gcosv,   &  ! cos/sin between model grid lines and NP direction at V point 
    36       gsinf, gcosf       ! cos/sin between model grid lines and NP direction at F point 
     29   PUBLIC   rot_rep   ! called in sbccpl, fldread, and cyclone 
     30   PUBLIC   geo2oce   ! called in sbccpl 
     31   PUBLIC   oce2geo   ! called in sbccpl 
     32   PUBLIC   obs_rot   ! called in obs_rot_vel and obs_write 
     33 
     34   !                                         ! cos/sin between model grid lines and NP direction 
     35   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsint, gcost   ! at T point 
     36   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinu, gcosu   ! at U point 
     37   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinv, gcosv   ! at V point 
     38   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinf, gcosf   ! at F point 
    3739 
    3840   LOGICAL ,              SAVE, DIMENSION(4)     ::   linit = .FALSE. 
    3941   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gsinlon, gcoslon, gsinlat, gcoslat 
    4042 
    41    LOGICAL ::   lmust_init = .TRUE.        !: used to initialize the cos/sin variables (se above) 
     43   LOGICAL ::   lmust_init = .TRUE.        !: used to initialize the cos/sin variables (see above) 
    4244 
    4345   !! * Substitutions 
     
    5052CONTAINS 
    5153 
    52    SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1,   & 
    53                        px2 , py2 ) 
    54       !!---------------------------------------------------------------------- 
    55       !!                  ***  ROUTINE repcmo  *** 
    56       !! 
    57       !! ** Purpose :   Change vector componantes from a geographic grid to a 
    58       !!      stretched coordinates grid. 
    59       !! 
    60       !! ** Method  :   Initialization of arrays at the first call. 
    61       !! 
    62       !! ** Action  : - px2 : first  componante (defined at u point) 
    63       !!              - py2 : second componante (defined at v point) 
    64       !!---------------------------------------------------------------------- 
    65       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pxu1, pyu1   ! geographic vector componantes at u-point 
    66       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pxv1, pyv1   ! geographic vector componantes at v-point 
    67       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2          ! i-componante (defined at u-point) 
    68       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   py2          ! j-componante (defined at v-point) 
    69       !!---------------------------------------------------------------------- 
    70        
    71       ! Change from geographic to stretched coordinate 
    72       ! ---------------------------------------------- 
    73       CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 ) 
    74       CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 ) 
    75        
    76    END SUBROUTINE repcmo 
    77  
    78  
    7954   SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot ) 
    8055      !!---------------------------------------------------------------------- 
     
    8358      !! ** Purpose :   Rotate the Repere: Change vector componantes between 
    8459      !!                geographic grid <--> stretched coordinates grid. 
    85       !! 
    86       !! History : 
    87       !!   9.2  !  07-04  (S. Masson)   
    88       !!                  (O. Marti ) Original code (repere and repcmo) 
    89       !!---------------------------------------------------------------------- 
    90       REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) ::   pxin, pyin   ! vector componantes 
    91       CHARACTER(len=1),             INTENT( IN ) ::   cd_type      ! define the nature of pt2d array grid-points 
    92       CHARACTER(len=5),             INTENT( IN ) ::   cdtodo       ! specify the work to do: 
    93       !!                                                           ! 'en->i' east-north componantes to model i componante 
    94       !!                                                           ! 'en->j' east-north componantes to model j componante 
    95       !!                                                           ! 'ij->e' model i-j componantes to east componante 
    96       !!                                                           ! 'ij->n' model i-j componantes to east componante 
    97       REAL(wp), DIMENSION(jpi,jpj), INTENT(out) ::   prot       
    98       !!---------------------------------------------------------------------- 
    99  
    100       ! Initialization of gsin* and gcos* at first call 
    101       ! ----------------------------------------------- 
    102  
    103       IF( lmust_init ) THEN 
     60      !!---------------------------------------------------------------------- 
     61      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pxin, pyin   ! vector componantes 
     62      CHARACTER(len=1),             INTENT(in   ) ::   cd_type      ! define the nature of pt2d array grid-points 
     63      CHARACTER(len=5),             INTENT(in   ) ::   cdtodo       ! type of transpormation: 
     64      !                                                             ! 'en->i' = east-north to i-component 
     65      !                                                             ! 'en->j' = east-north to j-component 
     66      !                                                             ! 'ij->e' = (i,j) components to east 
     67      !                                                             ! 'ij->n' = (i,j) components to north 
     68      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   prot       
     69      !!---------------------------------------------------------------------- 
     70      ! 
     71      IF( lmust_init ) THEN      ! at 1st call only: set  gsin. & gcos. 
    10472         IF(lwp) WRITE(numout,*) 
    105          IF(lwp) WRITE(numout,*) ' rot_rep : geographic <--> stretched' 
    106          IF(lwp) WRITE(numout,*) ' ~~~~~    coordinate transformation' 
     73         IF(lwp) WRITE(numout,*) ' rot_rep: coordinate transformation : geographic <==> model (i,j)-components' 
     74         IF(lwp) WRITE(numout,*) ' ~~~~~~~~    ' 
    10775         ! 
    10876         CALL angle       ! initialization of the transformation 
    10977         lmust_init = .FALSE. 
    11078      ENDIF 
    111        
    112       SELECT CASE (cdtodo) 
    113       CASE ('en->i')      ! 'en->i' est-north componantes to model i componante 
     79      ! 
     80      SELECT CASE( cdtodo )      ! type of rotation 
     81      ! 
     82      CASE( 'en->i' )                  ! east-north to i-component 
    11483         SELECT CASE (cd_type) 
    11584         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:) 
     
    11988         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 
    12089         END SELECT 
    121       CASE ('en->j')      ! 'en->j' est-north componantes to model j componante 
     90      CASE ('en->j')                   ! east-north to j-component 
    12291         SELECT CASE (cd_type) 
    12392         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:) 
     
    12796         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 
    12897         END SELECT 
    129       CASE ('ij->e')      ! 'ij->e' model i-j componantes to est componante 
     98      CASE ('ij->e')                   ! (i,j)-components to east 
    13099         SELECT CASE (cd_type) 
    131100         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:) 
     
    135104         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 
    136105         END SELECT 
    137       CASE ('ij->n')      ! 'ij->n' model i-j componantes to est componante 
     106      CASE ('ij->n')                   ! (i,j)-components to north  
    138107         SELECT CASE (cd_type) 
    139108         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:) 
     
    145114      CASE DEFAULT   ;   CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' ) 
    146115      END SELECT 
    147        
     116      ! 
    148117   END SUBROUTINE rot_rep 
    149118 
     
    155124      !! ** Purpose :   Compute angles between model grid lines and the North direction 
    156125      !! 
    157       !! ** Method  : 
    158       !! 
    159       !! ** Action  :   Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays: 
    160       !!      sinus and cosinus of the angle between the north-south axe and the  
    161       !!      j-direction at t, u, v and f-points 
    162       !! 
    163       !! History : 
    164       !!   7.0  !  96-07  (O. Marti )  Original code 
    165       !!   8.0  !  98-06  (G. Madec ) 
    166       !!   8.5  !  98-06  (G. Madec )  Free form, F90 + opt. 
    167       !!   9.2  !  07-04  (S. Masson)  Add T, F points and bugfix in cos lateral boundary 
    168       !!---------------------------------------------------------------------- 
    169       INTEGER ::   ji, jj   ! dummy loop indices 
    170       INTEGER ::   ierr     ! local integer 
    171       REAL(wp) ::   & 
    172          zlam, zphi,            &  ! temporary scalars 
    173          zlan, zphh,            &  !    "         " 
    174          zxnpt, zynpt, znnpt,   &  ! x,y components and norm of the vector: T point to North Pole 
    175          zxnpu, zynpu, znnpu,   &  ! x,y components and norm of the vector: U point to North Pole 
    176          zxnpv, zynpv, znnpv,   &  ! x,y components and norm of the vector: V point to North Pole 
    177          zxnpf, zynpf, znnpf,   &  ! x,y components and norm of the vector: F point to North Pole 
    178          zxvvt, zyvvt, znvvt,   &  ! x,y components and norm of the vector: between V points below and above a T point 
    179          zxffu, zyffu, znffu,   &  ! x,y components and norm of the vector: between F points below and above a U point 
    180          zxffv, zyffv, znffv,   &  ! x,y components and norm of the vector: between F points left  and right a V point 
    181          zxuuf, zyuuf, znuuf       ! x,y components and norm of the vector: between U points below and above a F point 
    182       !!---------------------------------------------------------------------- 
    183  
     126      !! ** Method  :   sinus and cosinus of the angle between the north-south axe  
     127      !!              and the j-direction at t, u, v and f-points 
     128      !!                dot and cross products are used to obtain cos and sin, resp. 
     129      !! 
     130      !! ** Action  : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf 
     131      !!---------------------------------------------------------------------- 
     132      INTEGER  ::   ji, jj   ! dummy loop indices 
     133      INTEGER  ::   ierr     ! local integer 
     134      REAL(wp) ::   zlam, zphi            ! local scalars 
     135      REAL(wp) ::   zlan, zphh            !   -      - 
     136      REAL(wp) ::   zxnpt, zynpt, znnpt   ! x,y components and norm of the vector: T point to North Pole 
     137      REAL(wp) ::   zxnpu, zynpu, znnpu   ! x,y components and norm of the vector: U point to North Pole 
     138      REAL(wp) ::   zxnpv, zynpv, znnpv   ! x,y components and norm of the vector: V point to North Pole 
     139      REAL(wp) ::   zxnpf, zynpf, znnpf   ! x,y components and norm of the vector: F point to North Pole 
     140      REAL(wp) ::   zxvvt, zyvvt, znvvt   ! x,y components and norm of the vector: between V points below and above a T point 
     141      REAL(wp) ::   zxffu, zyffu, znffu   ! x,y components and norm of the vector: between F points below and above a U point 
     142      REAL(wp) ::   zxffv, zyffv, znffv   ! x,y components and norm of the vector: between F points left  and right a V point 
     143      REAL(wp) ::   zxuuf, zyuuf, znuuf   ! x,y components and norm of the vector: between U points below and above a F point 
     144      !!---------------------------------------------------------------------- 
     145      ! 
    184146      ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj),   &  
    185147         &      gsinu(jpi,jpj), gcosu(jpi,jpj),   &  
     
    187149         &      gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr ) 
    188150      IF(lk_mpp)   CALL mpp_sum( ierr ) 
    189       IF( ierr /= 0 )   CALL ctl_stop('angle: unable to allocate arrays' ) 
    190  
     151      IF( ierr /= 0 )   CALL ctl_stop( 'angle: unable to allocate arrays' ) 
     152      ! 
    191153      ! ============================= ! 
    192154      ! Compute the cosinus and sinus ! 
    193155      ! ============================= ! 
    194156      ! (computation done on the north stereographic polar plane) 
    195  
     157      ! 
    196158      DO jj = 2, jpjm1 
    197159         DO ji = fs_2, jpi   ! vector opt. 
    198  
    199             ! north pole direction & modulous (at t-point) 
    200             zlam = glamt(ji,jj) 
     160            !                   
     161            zlam = glamt(ji,jj)     ! north pole direction & modulous (at t-point) 
    201162            zphi = gphit(ji,jj) 
    202163            zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    203164            zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    204165            znnpt = zxnpt*zxnpt + zynpt*zynpt 
    205  
    206             ! north pole direction & modulous (at u-point) 
    207             zlam = glamu(ji,jj) 
     166            ! 
     167            zlam = glamu(ji,jj)     ! north pole direction & modulous (at u-point) 
    208168            zphi = gphiu(ji,jj) 
    209169            zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    210170            zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    211171            znnpu = zxnpu*zxnpu + zynpu*zynpu 
    212  
    213             ! north pole direction & modulous (at v-point) 
    214             zlam = glamv(ji,jj) 
     172            ! 
     173            zlam = glamv(ji,jj)     ! north pole direction & modulous (at v-point) 
    215174            zphi = gphiv(ji,jj) 
    216175            zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    217176            zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    218177            znnpv = zxnpv*zxnpv + zynpv*zynpv 
    219  
    220             ! north pole direction & modulous (at f-point) 
    221             zlam = glamf(ji,jj) 
     178            ! 
     179            zlam = glamf(ji,jj)     ! north pole direction & modulous (at f-point) 
    222180            zphi = gphif(ji,jj) 
    223181            zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    224182            zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    225183            znnpf = zxnpf*zxnpf + zynpf*zynpf 
    226  
    227             ! j-direction: v-point segment direction (around t-point) 
    228             zlam = glamv(ji,jj  ) 
     184            ! 
     185            zlam = glamv(ji,jj  )   ! j-direction: v-point segment direction (around t-point) 
    229186            zphi = gphiv(ji,jj  ) 
    230187            zlan = glamv(ji,jj-1) 
     
    236193            znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  ) 
    237194            znvvt = MAX( znvvt, 1.e-14 ) 
    238  
    239             ! j-direction: f-point segment direction (around u-point) 
    240             zlam = glamf(ji,jj  ) 
     195            ! 
     196            zlam = glamf(ji,jj  )   ! j-direction: f-point segment direction (around u-point) 
    241197            zphi = gphif(ji,jj  ) 
    242198            zlan = glamf(ji,jj-1) 
     
    248204            znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu )  ) 
    249205            znffu = MAX( znffu, 1.e-14 ) 
    250  
    251             ! i-direction: f-point segment direction (around v-point) 
    252             zlam = glamf(ji  ,jj) 
     206            ! 
     207            zlam = glamf(ji  ,jj)   ! i-direction: f-point segment direction (around v-point) 
    253208            zphi = gphif(ji  ,jj) 
    254209            zlan = glamf(ji-1,jj) 
     
    260215            znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv )  ) 
    261216            znffv = MAX( znffv, 1.e-14 ) 
    262  
    263             ! j-direction: u-point segment direction (around f-point) 
    264             zlam = glamu(ji,jj+1) 
     217            ! 
     218            zlam = glamu(ji,jj+1)   ! j-direction: u-point segment direction (around f-point) 
    265219            zphi = gphiu(ji,jj+1) 
    266220            zlan = glamu(ji,jj  ) 
     
    272226            znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  ) 
    273227            znuuf = MAX( znuuf, 1.e-14 ) 
    274  
    275             ! cosinus and sinus using scalar and vectorial products 
     228            ! 
     229            !                       ! cosinus and sinus using dot and cross products 
    276230            gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt 
    277231            gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 
    278  
     232            ! 
    279233            gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu 
    280234            gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 
    281  
     235            ! 
    282236            gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf 
    283237            gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 
    284  
    285             ! (caution, rotation of 90 degres) 
     238            ! 
    286239            gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv 
    287             gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv 
    288  
     240            gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv     ! (caution, rotation of 90 degres) 
     241            ! 
    289242         END DO 
    290243      END DO 
     
    318271      ! Lateral boundary conditions ! 
    319272      ! =========================== ! 
    320  
    321       ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 
     273      !           ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 
    322274      CALL lbc_lnk( gcost, 'T', -1. )   ;   CALL lbc_lnk( gsint, 'T', -1. ) 
    323275      CALL lbc_lnk( gcosu, 'U', -1. )   ;   CALL lbc_lnk( gsinu, 'U', -1. ) 
    324276      CALL lbc_lnk( gcosv, 'V', -1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. ) 
    325277      CALL lbc_lnk( gcosf, 'F', -1. )   ;   CALL lbc_lnk( gsinf, 'F', -1. ) 
    326  
     278      ! 
    327279   END SUBROUTINE angle 
    328280 
    329281 
    330    SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid,     & 
    331                         pte, ptn ) 
     282   SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, pte, ptn ) 
    332283      !!---------------------------------------------------------------------- 
    333284      !!                    ***  ROUTINE geo2oce  *** 
     
    335286      !! ** Purpose : 
    336287      !! 
    337       !! ** Method  :   Change wind stress from geocentric to east/north 
    338       !! 
    339       !! History : 
    340       !!        !         (O. Marti)  Original code 
    341       !!        !  91-03  (G. Madec) 
    342       !!        !  92-07  (M. Imbard) 
    343       !!        !  99-11  (M. Imbard) NetCDF format with IOIPSL 
    344       !!        !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb 
    345       !!   8.5  !  02-06  (G. Madec)  F90: Free form 
    346       !!   3.0  !  07-08  (G. Madec)  geo2oce suppress lon/lat agruments 
     288      !! ** Method  :   Change a vector from geocentric to east/north  
     289      !! 
    347290      !!---------------------------------------------------------------------- 
    348291      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::  pxx, pyy, pzz 
    349292      CHARACTER(len=1)            , INTENT(in   ) ::  cgrid 
    350293      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::  pte, ptn 
    351       !! 
     294      ! 
    352295      REAL(wp), PARAMETER :: rpi = 3.141592653e0 
    353296      REAL(wp), PARAMETER :: rad = rpi / 180.e0 
     
    355298      INTEGER ::   ierr   ! local integer 
    356299      !!---------------------------------------------------------------------- 
    357  
     300      ! 
    358301      IF( .NOT. ALLOCATED( gsinlon ) ) THEN 
    359302         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   & 
     
    361304         IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    362305         IF( ierr /= 0 )   CALL ctl_stop('geo2oce: unable to allocate arrays' ) 
     306      ENDIF 
     307      ! 
     308      SELECT CASE( cgrid) 
     309      CASE ( 'T' )    
     310         ig = 1 
     311         IF( .NOT. linit(ig) ) THEN  
     312            gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 
     313            gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 
     314            gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 
     315            gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 
     316            linit(ig) = .TRUE. 
     317         ENDIF 
     318      CASE ( 'U' )    
     319         ig = 2 
     320         IF( .NOT. linit(ig) ) THEN  
     321            gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 
     322            gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 
     323            gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 
     324            gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 
     325            linit(ig) = .TRUE. 
     326         ENDIF 
     327      CASE ( 'V' )    
     328         ig = 3 
     329         IF( .NOT. linit(ig) ) THEN  
     330            gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 
     331            gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 
     332            gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 
     333            gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 
     334            linit(ig) = .TRUE. 
     335         ENDIF 
     336      CASE ( 'F' )    
     337         ig = 4 
     338         IF( .NOT. linit(ig) ) THEN  
     339            gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 
     340            gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 
     341            gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 
     342            gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 
     343            linit(ig) = .TRUE. 
     344         ENDIF 
     345      CASE default    
     346         WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 
     347         CALL ctl_stop( ctmp1 ) 
     348      END SELECT 
     349      ! 
     350      pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 
     351      ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx    & 
     352         &  - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy    & 
     353         &  + gcoslat(:,:,ig) * pzz 
     354      ! 
     355   END SUBROUTINE geo2oce 
     356 
     357 
     358   SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz ) 
     359      !!---------------------------------------------------------------------- 
     360      !!                    ***  ROUTINE oce2geo  *** 
     361      !!       
     362      !! ** Purpose : 
     363      !! 
     364      !! ** Method  :   Change vector from east/north to geocentric 
     365      !! 
     366      !! History :     ! (A. Caubel)  oce2geo - Original code 
     367      !!---------------------------------------------------------------------- 
     368      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN    ) ::  pte, ptn 
     369      CHARACTER(len=1)            , INTENT( IN    ) ::  cgrid 
     370      REAL(wp), DIMENSION(jpi,jpj), INTENT(   OUT ) ::  pxx , pyy , pzz 
     371      !! 
     372      REAL(wp), PARAMETER :: rpi = 3.141592653E0 
     373      REAL(wp), PARAMETER :: rad = rpi / 180.e0 
     374      INTEGER ::   ig     ! 
     375      INTEGER ::   ierr   ! local integer 
     376      !!---------------------------------------------------------------------- 
     377 
     378      IF( .NOT. ALLOCATED( gsinlon ) ) THEN 
     379         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   & 
     380            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 
     381         IF( lk_mpp    )   CALL mpp_sum( ierr ) 
     382         IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' ) 
    363383      ENDIF 
    364384 
     
    404424            CALL ctl_stop( ctmp1 ) 
    405425      END SELECT 
    406        
    407       pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 
    408       ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx    & 
    409             - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy    & 
    410             + gcoslat(:,:,ig) * pzz 
    411 !!$   ptv =   gcoslon(:,:,ig) * gcoslat(:,:,ig) * pxx    & 
    412 !!$         + gsinlon(:,:,ig) * gcoslat(:,:,ig) * pyy    & 
    413 !!$         + gsinlat(:,:,ig) * pzz 
    414       ! 
    415    END SUBROUTINE geo2oce 
    416  
    417    SUBROUTINE oce2geo ( pte, ptn, cgrid,     & 
    418                         pxx , pyy , pzz ) 
    419       !!---------------------------------------------------------------------- 
    420       !!                    ***  ROUTINE oce2geo  *** 
    421       !!       
    422       !! ** Purpose : 
    423       !! 
    424       !! ** Method  :   Change vector from east/north to geocentric 
    425       !! 
    426       !! History : 
    427       !!        !         (A. Caubel)  oce2geo - Original code 
    428       !!---------------------------------------------------------------------- 
    429       REAL(wp), DIMENSION(jpi,jpj), INTENT( IN    ) ::  pte, ptn 
    430       CHARACTER(len=1)            , INTENT( IN    ) ::  cgrid 
    431       REAL(wp), DIMENSION(jpi,jpj), INTENT(   OUT ) ::  pxx , pyy , pzz 
    432       !! 
    433       REAL(wp), PARAMETER :: rpi = 3.141592653E0 
    434       REAL(wp), PARAMETER :: rad = rpi / 180.e0 
    435       INTEGER ::   ig     ! 
    436       INTEGER ::   ierr   ! local integer 
    437       !!---------------------------------------------------------------------- 
    438  
    439       IF( .NOT. ALLOCATED( gsinlon ) ) THEN 
    440          ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   & 
    441             &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 
    442          IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    443          IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' ) 
    444       ENDIF 
    445  
    446       SELECT CASE( cgrid) 
    447          CASE ( 'T' )    
    448             ig = 1 
    449             IF( .NOT. linit(ig) ) THEN  
    450                gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 
    451                gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 
    452                gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 
    453                gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 
    454                linit(ig) = .TRUE. 
    455             ENDIF 
    456          CASE ( 'U' )    
    457             ig = 2 
    458             IF( .NOT. linit(ig) ) THEN  
    459                gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 
    460                gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 
    461                gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 
    462                gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 
    463                linit(ig) = .TRUE. 
    464             ENDIF 
    465          CASE ( 'V' )    
    466             ig = 3 
    467             IF( .NOT. linit(ig) ) THEN  
    468                gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 
    469                gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 
    470                gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 
    471                gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 
    472                linit(ig) = .TRUE. 
    473             ENDIF 
    474          CASE ( 'F' )    
    475             ig = 4 
    476             IF( .NOT. linit(ig) ) THEN  
    477                gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 
    478                gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 
    479                gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 
    480                gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 
    481                linit(ig) = .TRUE. 
    482             ENDIF 
    483          CASE default    
    484             WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 
    485             CALL ctl_stop( ctmp1 ) 
    486       END SELECT 
    487  
    488        pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn  
    489        pyy =   gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn 
    490        pzz =   gcoslat(:,:,ig) * ptn 
    491  
    492        
     426      ! 
     427      pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn  
     428      pyy =   gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn 
     429      pzz =   gcoslat(:,:,ig) * ptn 
     430      ! 
    493431   END SUBROUTINE oce2geo 
    494432 
    495433 
    496    SUBROUTINE repere ( px1, py1, px2, py2, kchoix, cd_type ) 
    497       !!---------------------------------------------------------------------- 
    498       !!                 ***  ROUTINE repere  *** 
    499       !!         
    500       !! ** Purpose :   Change vector componantes between a geopgraphic grid  
    501       !!      and a stretched coordinates grid. 
    502       !! 
    503       !! ** Method  :    
    504       !! 
    505       !! ** Action  : 
    506       !! 
    507       !! History : 
    508       !!        !  89-03  (O. Marti)  original code 
    509       !!        !  92-02  (M. Imbard) 
    510       !!        !  93-03  (M. Guyon)  symetrical conditions 
    511       !!        !  98-05  (B. Blanke) 
    512       !!   8.5  !  02-08  (G. Madec)  F90: Free form 
    513       !!---------------------------------------------------------------------- 
    514       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   px1, py1   ! two horizontal components to be rotated 
    515       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2, py2   ! the two horizontal components in the model repere 
    516       INTEGER , INTENT(in   )                     ::   kchoix     ! type of transformation 
    517       !                                                           ! = 1 change from geographic to model grid. 
    518       !                                                           ! =-1 change from model to geographic grid 
    519       CHARACTER(len=1), INTENT(in   ), OPTIONAL   ::   cd_type    ! define the nature of pt2d array grid-points 
    520       ! 
    521       CHARACTER(len=1) ::   cl_type      ! define the nature of pt2d array grid-points (T point by default) 
    522       !!---------------------------------------------------------------------- 
    523  
    524       cl_type = 'T' 
    525       IF( PRESENT(cd_type) )   cl_type = cd_type 
    526          ! 
    527       SELECT CASE (kchoix) 
    528       CASE ( 1)      ! change from geographic to model grid. 
    529          CALL rot_rep( px1, py1, cl_type, 'en->i', px2 ) 
    530          CALL rot_rep( px1, py1, cl_type, 'en->j', py2 ) 
    531       CASE (-1)      ! change from model to geographic grid 
    532          CALL rot_rep( px1, py1, cl_type, 'ij->e', px2 ) 
    533          CALL rot_rep( px1, py1, cl_type, 'ij->n', py2 ) 
    534       CASE DEFAULT   ;   CALL ctl_stop( 'repere: Syntax Error in the definition of kchoix (1 OR -1' ) 
    535       END SELECT 
    536        
    537    END SUBROUTINE repere 
    538  
    539  
    540    SUBROUTINE obs_rot ( psinu, pcosu, psinv, pcosv ) 
     434   SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv ) 
    541435      !!---------------------------------------------------------------------- 
    542436      !!                  ***  ROUTINE obs_rot  *** 
     
    546440      !!                current at observation points 
    547441      !! 
    548       !! History : 
    549       !!   9.2  !  09-02  (K. Mogensen) 
     442      !! History :  9.2  !  09-02  (K. Mogensen) 
    550443      !!---------------------------------------------------------------------- 
    551444      REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT )::   psinu, pcosu, psinv, pcosv   ! copy of data 
    552445      !!---------------------------------------------------------------------- 
    553  
     446      ! 
    554447      ! Initialization of gsin* and gcos* at first call 
    555448      ! ----------------------------------------------- 
    556  
    557449      IF( lmust_init ) THEN 
    558450         IF(lwp) WRITE(numout,*) 
    559451         IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched' 
    560452         IF(lwp) WRITE(numout,*) ' ~~~~~~~   coordinate transformation' 
    561  
    562453         CALL angle       ! initialization of the transformation 
    563454         lmust_init = .FALSE. 
    564  
    565455      ENDIF 
    566  
     456      ! 
    567457      psinu(:,:) = gsinu(:,:) 
    568458      pcosu(:,:) = gcosu(:,:) 
    569459      psinv(:,:) = gsinv(:,:) 
    570460      pcosv(:,:) = gcosv(:,:) 
    571  
     461      ! 
    572462   END SUBROUTINE obs_rot 
    573463 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r5845 r5883  
    1717 
    1818   !!---------------------------------------------------------------------- 
    19    !!   sbc_init       : read namsbc namelist 
    20    !!   sbc            : surface ocean momentum, heat and freshwater boundary conditions 
     19   !!   sbc_init      : read namsbc namelist 
     20   !!   sbc           : surface ocean momentum, heat and freshwater boundary conditions 
    2121   !!---------------------------------------------------------------------- 
    22    USE oce              ! ocean dynamics and tracers 
    23    USE dom_oce          ! ocean space and time domain 
    24    USE phycst           ! physical constants 
    25    USE sbc_oce          ! Surface boundary condition: ocean fields 
    26    USE trc_oce          ! shared ocean-passive tracers variables 
    27    USE sbc_ice          ! Surface boundary condition: ice fields 
    28    USE sbcdcy           ! surface boundary condition: diurnal cycle 
    29    USE sbcssm           ! surface boundary condition: sea-surface mean variables 
    30    USE sbcana           ! surface boundary condition: analytical formulation 
    31    USE sbcflx           ! surface boundary condition: flux formulation 
    32    USE sbcblk_clio      ! surface boundary condition: bulk formulation : CLIO 
    33    USE sbcblk_core      ! surface boundary condition: bulk formulation : CORE 
    34    USE sbcblk_mfs       ! surface boundary condition: bulk formulation : MFS 
    35    USE sbcice_if        ! surface boundary condition: ice-if sea-ice model 
    36    USE sbcice_lim       ! surface boundary condition: LIM 3.0 sea-ice model 
    37    USE sbcice_lim_2     ! surface boundary condition: LIM 2.0 sea-ice model 
    38    USE sbcice_cice      ! surface boundary condition: CICE    sea-ice model 
    39    USE sbccpl           ! surface boundary condition: coupled florulation 
    40    USE cpl_oasis3       ! OASIS routines for coupling 
    41    USE sbcssr           ! surface boundary condition: sea surface restoring 
    42    USE sbcrnf           ! surface boundary condition: runoffs 
    43    USE sbcisf           ! surface boundary condition: ice shelf 
    44    USE sbcfwb           ! surface boundary condition: freshwater budget 
    45    USE closea           ! closed sea 
    46    USE icbstp           ! Icebergs! 
    47  
    48    USE prtctl           ! Print control                    (prt_ctl routine) 
    49    USE iom              ! IOM library 
    50    USE in_out_manager   ! I/O manager 
    51    USE lib_mpp          ! MPP library 
    52    USE timing           ! Timing 
    53    USE sbcwave          ! Wave module 
    54    USE bdy_par          ! Require lk_bdy 
     22   USE oce            ! ocean dynamics and tracers 
     23   USE dom_oce        ! ocean space and time domain 
     24   USE phycst         ! physical constants 
     25   USE sbc_oce        ! Surface boundary condition: ocean fields 
     26   USE trc_oce        ! shared ocean-passive tracers variables 
     27   USE sbc_ice        ! Surface boundary condition: ice fields 
     28   USE sbcdcy         ! surface boundary condition: diurnal cycle 
     29   USE sbcssm         ! surface boundary condition: sea-surface mean variables 
     30   USE sbcana         ! surface boundary condition: analytical formulation 
     31   USE sbcflx         ! surface boundary condition: flux formulation 
     32   USE sbcblk_clio    ! surface boundary condition: bulk formulation : CLIO 
     33   USE sbcblk_core    ! surface boundary condition: bulk formulation : CORE 
     34   USE sbcblk_mfs     ! surface boundary condition: bulk formulation : MFS 
     35   USE sbcice_if      ! surface boundary condition: ice-if sea-ice model 
     36   USE sbcice_lim     ! surface boundary condition: LIM 3.0 sea-ice model 
     37   USE sbcice_lim_2   ! surface boundary condition: LIM 2.0 sea-ice model 
     38   USE sbcice_cice    ! surface boundary condition: CICE    sea-ice model 
     39   USE sbccpl         ! surface boundary condition: coupled florulation 
     40   USE cpl_oasis3     ! OASIS routines for coupling 
     41   USE sbcssr         ! surface boundary condition: sea surface restoring 
     42   USE sbcrnf         ! surface boundary condition: runoffs 
     43   USE sbcisf         ! surface boundary condition: ice shelf 
     44   USE sbcfwb         ! surface boundary condition: freshwater budget 
     45   USE closea         ! closed sea 
     46   USE icbstp         ! Icebergs 
     47   USE traqsr         ! active tracers: light penetration 
     48   USE sbcwave        ! Wave module 
     49   USE bdy_par        ! Require lk_bdy 
     50   ! 
     51   USE prtctl         ! Print control                    (prt_ctl routine) 
     52   USE iom            ! IOM library 
     53   USE in_out_manager ! I/O manager 
     54   USE lib_mpp        ! MPP library 
     55   USE timing         ! Timing 
    5556 
    5657   IMPLICIT NONE 
     
    8384      INTEGER ::   icpt   ! local integer 
    8485      !! 
    85       NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   & 
    86          &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   & 
    87          &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   & 
    88          &             nn_lsm    , nn_limflx , nn_components, ln_cpl 
     86      NAMELIST/namsbc/ nn_fsbc  , ln_ana   , ln_flx, ln_blk_clio, ln_blk_core, ln_blk_mfs,   & 
     87         &             ln_cpl   , ln_mixcpl, nn_components      , nn_limflx  ,               & 
     88         &             ln_traqsr, ln_dm2dc ,                                                 &   
     89         &             nn_ice   , nn_ice_embd,                                               & 
     90         &             ln_rnf   , ln_ssr   , nn_isf   , nn_fwb    , ln_apr_dyn,              & 
     91         &             ln_wave  ,                                                            & 
     92         &             nn_lsm    
    8993      INTEGER  ::   ios 
    9094      INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm 
    9195      LOGICAL  ::   ll_purecpl 
    9296      !!---------------------------------------------------------------------- 
    93  
     97      ! 
    9498      IF(lwp) THEN 
    9599         WRITE(numout,*) 
     
    97101         WRITE(numout,*) '~~~~~~~~ ' 
    98102      ENDIF 
    99  
     103      ! 
    100104      REWIND( numnam_ref )              ! Namelist namsbc in reference namelist : Surface boundary 
    101105      READ  ( numnam_ref, namsbc, IOSTAT = ios, ERR = 901) 
    102 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp ) 
    103  
     106901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp ) 
     107      ! 
    104108      REWIND( numnam_cfg )              ! Namelist namsbc in configuration namelist : Parameters of the run 
    105109      READ  ( numnam_cfg, namsbc, IOSTAT = ios, ERR = 902 ) 
    106 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp ) 
    107       IF(lwm) WRITE ( numond, namsbc ) 
    108  
     110902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp ) 
     111      IF(lwm) WRITE( numond, namsbc ) 
     112      ! 
    109113      !                          ! overwrite namelist parameter using CPP key information 
    110114      IF( Agrif_Root() ) THEN                ! AGRIF zoom 
     
    117121          nn_ice      =   0 
    118122      ENDIF 
    119  
     123      ! 
    120124      IF(lwp) THEN               ! Control print 
    121125         WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)' 
    122126         WRITE(numout,*) '           frequency update of sbc (and ice)             nn_fsbc     = ', nn_fsbc 
    123          WRITE(numout,*) '           Type of sbc : ' 
    124          WRITE(numout,*) '              analytical formulation                     ln_ana      = ', ln_ana 
    125          WRITE(numout,*) '              flux       formulation                     ln_flx      = ', ln_flx 
    126          WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_clio = ', ln_blk_clio 
    127          WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core = ', ln_blk_core 
    128          WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs 
    129          WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl 
    130          WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl   = ', ln_mixcpl 
    131          WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis 
    132          WRITE(numout,*) '              components of your executable            nn_components = ', nn_components 
    133          WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx 
     127         WRITE(numout,*) '           Type of air-sea fluxes : ' 
     128         WRITE(numout,*) '              analytical formulation                     ln_ana        = ', ln_ana 
     129         WRITE(numout,*) '              flux       formulation                     ln_flx        = ', ln_flx 
     130         WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_clio   = ', ln_blk_clio 
     131         WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core   = ', ln_blk_core 
     132         WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs    = ', ln_blk_mfs 
     133         WRITE(numout,*) '           Type of coupling (Ocean/Ice/Atmosphere) : ' 
     134         WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl 
     135         WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl     = ', ln_mixcpl 
     136         WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis      = ', lk_oasis 
     137         WRITE(numout,*) '              components of your executable              nn_components = ', nn_components 
     138         WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx     = ', nn_limflx 
     139         WRITE(numout,*) '           Sea-ice : ' 
     140         WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice        = ', nn_ice  
     141         WRITE(numout,*) '              ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd   = ', nn_ice_embd 
    134142         WRITE(numout,*) '           Misc. options of sbc : ' 
    135          WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn 
    136          WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice      = ', nn_ice  
    137          WRITE(numout,*) '              ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd = ', nn_ice_embd 
    138          WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc  
    139          WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf 
    140          WRITE(numout,*) '              iceshelf formulation                       nn_isf      = ', nn_isf 
    141          WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr      = ', ln_ssr 
    142          WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb      = ', nn_fwb 
    143          WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea 
    144          WRITE(numout,*) '              n. of iterations if land-sea-mask applied  nn_lsm      = ', nn_lsm 
    145       ENDIF 
    146  
    147       ! LIM3 Multi-category heat flux formulation 
    148       SELECT CASE ( nn_limflx) 
    149       CASE ( -1 ) 
    150          IF(lwp) WRITE(numout,*) '              Use of per-category fluxes (nn_limflx = -1) ' 
    151       CASE ( 0  ) 
    152          IF(lwp) WRITE(numout,*) '              Average per-category fluxes (nn_limflx = 0) '  
    153       CASE ( 1  ) 
    154          IF(lwp) WRITE(numout,*) '              Average then redistribute per-category fluxes (nn_limflx = 1) ' 
    155       CASE ( 2  ) 
    156          IF(lwp) WRITE(numout,*) '              Redistribute a single flux over categories (nn_limflx = 2) ' 
    157       END SELECT 
    158       ! 
    159       IF ( nn_components /= jp_iam_nemo .AND. .NOT. lk_oasis )   & 
    160          &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' ) 
    161       IF ( nn_components == jp_iam_opa .AND. ln_cpl )   & 
    162          &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' ) 
    163       IF ( nn_components == jp_iam_opa .AND. ln_mixcpl )   & 
    164          &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' ) 
    165       IF ( ln_cpl .AND. .NOT. lk_oasis )    & 
    166          &      CALL ctl_stop( 'STOP', 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' ) 
     143         WRITE(numout,*) '              Light penetration in temperature Eq.       ln_traqsr     = ', ln_traqsr 
     144         WRITE(numout,*) '                 daily mean to diurnal cycle qsr            ln_dm2dc   = ', ln_dm2dc  
     145         WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr        = ', ln_ssr 
     146         WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb        = ', nn_fwb 
     147         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn    = ', ln_apr_dyn 
     148         WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf        = ', ln_rnf 
     149         WRITE(numout,*) '              iceshelf formulation                       nn_isf        = ', nn_isf 
     150         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea     = ', nn_closea 
     151         WRITE(numout,*) '              nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm 
     152         WRITE(numout,*) '              surface wave                               ln_wave       = ', ln_wave   
     153      ENDIF 
     154      ! 
     155      IF(lwp) THEN 
     156         WRITE(numout,*) 
     157         SELECT CASE ( nn_limflx )        ! LIM3 Multi-category heat flux formulation 
     158         CASE ( -1 )   ;   WRITE(numout,*) '   LIM3: use per-category fluxes (nn_limflx = -1) ' 
     159         CASE ( 0  )   ;   WRITE(numout,*) '   LIM3: use average per-category fluxes (nn_limflx = 0) '  
     160         CASE ( 1  )   ;   WRITE(numout,*) '   LIM3: use average then redistribute per-category fluxes (nn_limflx = 1) ' 
     161         CASE ( 2  )   ;   WRITE(numout,*) '   LIM3: Redistribute a single flux over categories (nn_limflx = 2) ' 
     162         END SELECT 
     163      ENDIF 
     164      ! 
     165      IF( nn_components /= jp_iam_nemo .AND. .NOT. lk_oasis )   & 
     166         &      CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' ) 
     167      IF( nn_components == jp_iam_opa .AND. ln_cpl )   & 
     168         &      CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' ) 
     169      IF( nn_components == jp_iam_opa .AND. ln_mixcpl )   & 
     170         &      CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' ) 
     171      IF( ln_cpl .AND. .NOT. lk_oasis )    & 
     172         &      CALL ctl_stop( 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' ) 
    167173      IF( ln_mixcpl .AND. .NOT. lk_oasis )    & 
    168174         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires the cpp key key_oasis3' ) 
     
    176182 
    177183      !                          ! Checks: 
    178       IF( nn_isf .EQ. 0 ) THEN                      ! variable initialisation if no ice shelf  
     184      IF( nn_isf == 0 ) THEN                       ! variable initialisation if no ice shelf  
    179185         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_isf arrays' ) 
    180          fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp 
    181          risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp 
    182          rdivisf       = 0.0_wp 
     186         fwfisf  (:,:)   = 0._wp   ;   fwfisf_b  (:,:)   = 0._wp 
     187         risf_tsc(:,:,:) = 0._wp   ;   risf_tsc_b(:,:,:) = 0._wp 
     188         rdivisf       = 0._wp 
    183189      END IF 
    184       IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero 
    185  
    186       sfx(:,:) = 0.0_wp                            ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero)  
     190      IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! no ice in the domain, ice fraction is always zero 
     191 
     192      sfx(:,:) = 0._wp                             ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero)  
    187193                                                   ! only if sea-ice is present 
    188194  
    189       fmmflx(:,:) = 0.0_wp                        ! freezing-melting array initialisation 
     195      fmmflx(:,:) = 0._wp                          ! freezing-melting array initialisation 
    190196       
    191       taum(:,:) = 0.0_wp                           ! Initialise taum for use in gls in case of reduced restart 
     197      taum(:,:) = 0._wp                            ! Initialise taum for use in gls in case of reduced restart 
    192198 
    193199      !                                            ! restartability    
     
    212218         &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' ) 
    213219       
    214       IF ( ln_wave ) THEN 
    215       !Activated wave module but neither drag nor stokes drift activated 
    216          IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN 
    217             CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' ) 
    218       !drag coefficient read from wave model definable only with mfs bulk formulae and core  
    219          ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN        
    220              CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 
    221          ENDIF 
    222       ELSE 
    223       IF ( ln_cdgw .OR. ln_sdw  )                                                           &  
    224          &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
    225          &                  'with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ') 
    226       ENDIF  
    227220      !                          ! Choice of the Surface Boudary Condition (set nsbc) 
    228221      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl 
     
    243236      IF(lwp) THEN 
    244237         WRITE(numout,*) 
    245          IF( nsbc == jp_gyre    )   WRITE(numout,*) '              GYRE analytical formulation' 
    246          IF( nsbc == jp_ana     )   WRITE(numout,*) '              analytical formulation' 
    247          IF( nsbc == jp_flx     )   WRITE(numout,*) '              flux formulation' 
    248          IF( nsbc == jp_clio    )   WRITE(numout,*) '              CLIO bulk formulation' 
    249          IF( nsbc == jp_core    )   WRITE(numout,*) '              CORE bulk formulation' 
    250          IF( nsbc == jp_purecpl )   WRITE(numout,*) '              pure coupled formulation' 
    251          IF( nsbc == jp_mfs     )   WRITE(numout,*) '              MFS Bulk formulation' 
    252          IF( nsbc == jp_none    )   WRITE(numout,*) '              OPA coupled to SAS via oasis' 
    253          IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed formulation' 
     238         SELECT CASE( nsbc ) 
     239         CASE( jp_gyre    )   ;   WRITE(numout,*) '   GYRE analytical formulation' 
     240         CASE( jp_ana     )   ;   WRITE(numout,*) '   analytical formulation' 
     241         CASE( jp_flx     )   ;   WRITE(numout,*) '   flux formulation' 
     242         CASE( jp_clio    )   ;   WRITE(numout,*) '   CLIO bulk formulation' 
     243         CASE( jp_core    )   ;   WRITE(numout,*) '   CORE bulk formulation' 
     244         CASE( jp_purecpl )   ;   WRITE(numout,*) '   pure coupled formulation' 
     245         CASE( jp_mfs     )   ;   WRITE(numout,*) '   MFS Bulk formulation' 
     246         CASE( jp_none    )   ;   WRITE(numout,*) '   OPA coupled to SAS via oasis' 
     247            IF( ln_mixcpl )       WRITE(numout,*) '       + forced-coupled mixed formulation' 
     248         END SELECT 
    254249         IF( nn_components/= jp_iam_nemo )  & 
    255             &                       WRITE(numout,*) '              + OASIS coupled SAS' 
     250            &                     WRITE(numout,*) '       + OASIS coupled SAS' 
    256251      ENDIF 
    257252      ! 
    258253      IF( lk_oasis )   CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step 
    259254      !                                             !                                            (2) the use of nn_fsbc 
    260  
    261 !     nn_fsbc initialization if OPA-SAS coupling via OASIS 
    262 !     sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly 
    263       IF ( nn_components /= jp_iam_nemo ) THEN 
    264          IF ( nn_components == jp_iam_opa ) nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt) 
    265          IF ( nn_components == jp_iam_sas ) nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt) 
     255      !     nn_fsbc initialization if OPA-SAS coupling via OASIS 
     256      !     sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly 
     257      IF( nn_components /= jp_iam_nemo ) THEN 
     258         IF( nn_components == jp_iam_opa )   nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt) 
     259         IF( nn_components == jp_iam_sas )   nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt) 
    266260         ! 
    267261         IF(lwp)THEN 
     
    271265         ENDIF 
    272266      ENDIF 
    273  
     267      ! 
    274268      IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   & 
    275269          MOD( nstock             , nn_fsbc) /= 0 ) THEN  
     
    284278      IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) )  < 8 ) )   & 
    285279         &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' ) 
    286  
    287                                CALL sbc_ssm_init               ! Sea-surface mean fields initialisation 
    288       ! 
    289       IF( ln_ssr           )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation 
    290       ! 
    291                                CALL sbc_rnf_init               ! Runof initialisation 
    292       ! 
    293       IF( nn_ice == 3      )   CALL sbc_lim_init               ! LIM3 initialisation 
    294  
    295       IF( nn_ice == 4      )   CALL cice_sbc_init( nsbc )      ! CICE initialisation 
    296        
     280      ! 
     281                          CALL sbc_ssm_init               ! Sea-surface mean fields initialisation 
     282      ! 
     283      IF( ln_ssr      )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation 
     284      ! 
     285                          CALL sbc_rnf_init               ! Runof initialisation 
     286      ! 
     287      IF( nn_ice == 3 )   CALL sbc_lim_init               ! LIM3 initialisation 
     288      ! 
     289      IF( nn_ice == 4 )   CALL cice_sbc_init( nsbc )      ! CICE initialisation 
     290      ! 
    297291   END SUBROUTINE sbc_init 
    298292 
     
    325319         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields 
    326320         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine) 
    327          ! The 3D heat content due to qsr forcing is treated in traqsr 
    328          ! qsr_b (:,:) = qsr (:,:) 
    329          emp_b(:,:) = emp(:,:) 
    330          sfx_b(:,:) = sfx(:,:) 
     321         emp_b (:,:) = emp (:,:) 
     322         sfx_b (:,:) = sfx (:,:) 
    331323      ENDIF 
    332324      !                                            ! ---------------------------------------- ! 
     
    334326      !                                            ! ---------------------------------------- ! 
    335327      ! 
    336       IF( nn_components /= jp_iam_sas )   CALL sbc_ssm( kt )   ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 
     328      IF( nn_components /= jp_iam_sas )   CALL sbc_ssm ( kt )  ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 
    337329      !                                                        ! averaged over nf_sbc time-step 
    338  
    339       IF (ln_wave) CALL sbc_wave( kt ) 
     330      IF( ln_wave                     )   CALL sbc_wave( kt )  ! surface waves 
     331       
     332       
    340333                                                   !==  sbc formulation  ==! 
    341334                                                             
     
    355348      CASE( jp_mfs   )   ;   CALL sbc_blk_mfs ( kt )                    ! bulk formulation : MFS for the ocean 
    356349      CASE( jp_none  )  
    357          IF( nn_components == jp_iam_opa ) & 
    358                              CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS 
     350         IF( nn_components == jp_iam_opa )   & 
     351            &                CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS 
    359352      END SELECT 
    360353 
    361354      IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing 
    362355 
    363  
     356      ! 
    364357      !                                            !==  Misc. Options  ==! 
    365        
     358      ! 
    366359      SELECT CASE( nn_ice )                                       ! Update heat and freshwater fluxes over sea-ice areas 
    367360      CASE(  1 )   ;         CALL sbc_ice_if   ( kt )                ! Ice-cover climatology ("Ice-if" model) 
     
    428421         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx  ) 
    429422      ENDIF 
    430  
    431423      !                                                ! ---------------------------------------- ! 
    432424      !                                                !        Outputs and control print         ! 
     
    450442      ! 
    451443      IF(ln_ctl) THEN         ! print mean trends (used for debugging) 
    452          CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 ) 
    453          CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 ) 
    454          CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf  - : ', mask1=tmask, ovlap=1 ) 
     444         CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 ) 
     445         CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf - : ', mask1=tmask, ovlap=1 ) 
     446         CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf - : ', mask1=tmask, ovlap=1 ) 
    455447         CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 ) 
    456448         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 ) 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r5866 r5883  
    125125      IF(   ln_rnf_sal   )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required 
    126126      ! 
     127!!gm BUG :  TO BE REMOVED !! 
    127128      ! Runoff reduction only associated to the ORCA2_LIM configuration 
    128129      ! when reading the NetCDF file runoff_1m_nomask.nc 
     
    132133         END WHERE 
    133134      ENDIF 
     135!!gm end 
    134136      ! 
    135137      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r5845 r5883  
    7171      REAL(wp), DIMENSION(:,:,:), POINTER ::   zusd_t, zvsd_t, ze3hdiv   ! 3D workspace 
    7272      !! 
    73       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn 
     73      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn, ln_cdgw , ln_sdw 
    7474      !!--------------------------------------------------------------------- 
    7575      ! 
     
    8080         READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
    8181901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
    82  
     82         ! 
    8383         REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
    8484         READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
     
    8686         IF(lwm) WRITE ( numond, namsbc_wave ) 
    8787         ! 
    88          IF ( ln_cdgw ) THEN 
     88         IF(lwp) THEN               ! Control print 
     89            WRITE(numout,*) '        Namelist namsbc_wave : surface wave setting'  
     90            WRITE(numout,*) '           wave drag coefficient                      ln_cdgw  = ', ln_cdgw   
     91            WRITE(numout,*) '           wave stokes drift                          ln_sdw   = ', ln_sdw 
     92         ENDIF 
     93         ! 
     94         IF( .NOT.( ln_cdgw .OR. ln_sdw ) )    & 
     95            &  CALL ctl_warn( 'ln_sbcwave=T but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' ) 
     96         IF( ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )   &        
     97            &  CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 
     98         ! 
     99         IF( ln_cdgw ) THEN 
    89100            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    90101            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     
    96107            cdn_wave(:,:) = 0.0 
    97108         ENDIF 
    98          IF ( ln_sdw ) THEN 
     109         IF( ln_sdw ) THEN 
    99110            slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn 
    100111            ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
     
    114125      ENDIF 
    115126      ! 
    116       IF ( ln_cdgw ) THEN              !==  Neutral drag coefficient  ==! 
     127      IF( ln_cdgw ) THEN               !==  Neutral drag coefficient  ==! 
    117128         CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing 
    118129         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    119130      ENDIF 
    120131      ! 
    121       IF ( ln_sdw )  THEN              !==  Computation of the 3d Stokes Drift  ==! 
     132      IF( ln_sdw )  THEN               !==  Computation of the 3d Stokes Drift  ==! 
     133         ! 
     134         CALL wrk_alloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv ) 
    122135         ! 
    123136         CALL fld_read( kt, nn_fsbc, sf_sd )    !* read drag coefficient from external forcing 
    124137         ! 
    125          ! 
    126          CALL wrk_alloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv ) 
    127          !                                      !* distribute it on the vertical 
    128          DO jk = 1, jpkm1 
     138         DO jk = 1, jpkm1                       !* distribute it on the vertical 
    129139            zusd_t(:,:,jk) = sf_sd(jp_usd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * gdept_n(:,:,jk) ) 
    130140            zvsd_t(:,:,jk) = sf_sd(jp_vsd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * gdept_n(:,:,jk) ) 
    131141         END DO 
    132          !                                      !* interpolate the stokes drift from t-point to u- and v-points 
    133          DO jk = 1, jpkm1 
     142         DO jk = 1, jpkm1                       !* interpolate the stokes drift from t-point to u- and v-points 
    134143            DO jj = 1, jpjm1 
    135144               DO ji = 1, jpim1 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r5866 r5883  
    142142      CASE ( np_QCK )                                    ! QUICKEST 
    143143         CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts                     ) 
    144       ! 
    145144      END SELECT 
    146145      ! 
    147       !                                              ! print mean trends (used for debugging) 
     146      !                                         ! print mean trends (used for debugging) 
    148147      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv  - Ta: ', mask1=tmask,               & 
    149148         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     
    173172      ! 
    174173      !                                !==  Namelist  ==! 
    175       ! 
    176174      REWIND( numnam_ref )                   ! Namelist namtra_adv in reference namelist : Tracer advection scheme 
    177175      READ  ( numnam_ref, namtra_adv, IOSTAT = ios, ERR = 901) 
    178 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in reference namelist', lwp ) 
     176901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_adv in reference namelist', lwp ) 
    179177      ! 
    180178      REWIND( numnam_cfg )                   ! Namelist namtra_adv in configuration namelist : Tracer advection scheme 
    181179      READ  ( numnam_cfg, namtra_adv, IOSTAT = ios, ERR = 902 ) 
    182 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in configuration namelist', lwp ) 
    183       IF(lwm) WRITE ( numond, namtra_adv ) 
    184  
     180902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_adv in configuration namelist', lwp ) 
     181      IF(lwm) WRITE( numond, namtra_adv ) 
     182      ! 
    185183      IF(lwp) THEN                           ! Namelist print 
    186184         WRITE(numout,*) 
     
    201199         WRITE(numout,*) '      QUICKEST scheme                           ln_traadv_qck = ', ln_traadv_qck 
    202200      ENDIF 
    203  
     201      ! 
    204202      ioptio = 0                       !==  Parameter control  ==! 
    205203      IF( ln_traadv_cen )   ioptio = ioptio + 1 
     
    252250      IF( ln_traadv_ubs                      )   nadv = np_UBS 
    253251      IF( ln_traadv_qck                      )   nadv = np_QCK 
    254  
     252      ! 
    255253      IF(lwp) THEN                           ! Print the choice 
    256254         WRITE(numout,*) 
    257          IF( nadv == np_NO_adv  )   WRITE(numout,*) '         NO T-S advection' 
    258          IF( nadv == np_CEN     )   WRITE(numout,*) '         CEN      scheme is used. Horizontal order: ', nn_cen_h,   & 
    259             &                                                                        ' Vertical   order: ', nn_cen_v 
    260          IF( nadv == np_FCT     )   WRITE(numout,*) '         FCT      scheme is used. Horizontal order: ', nn_fct_h,   & 
    261             &                                                                        ' Vertical   order: ', nn_fct_v 
    262          IF( nadv == np_FCT_zts )   WRITE(numout,*) '         use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping' 
    263          IF( nadv == np_MUS     )   WRITE(numout,*) '         MUSCL    scheme is used' 
    264          IF( nadv == np_UBS     )   WRITE(numout,*) '         UBS      scheme is used' 
    265          IF( nadv == np_QCK     )   WRITE(numout,*) '         QUICKEST scheme is used' 
     255         SELECT CASE ( nadv ) 
     256         CASE( np_NO_adv  )   ;   WRITE(numout,*) '         NO T-S advection' 
     257         CASE( np_CEN     )   ;   WRITE(numout,*) '         CEN      scheme is used. Horizontal order: ', nn_cen_h,   & 
     258            &                                                                     ' Vertical   order: ', nn_cen_v 
     259         CASE( np_FCT     )   ;   WRITE(numout,*) '         FCT      scheme is used. Horizontal order: ', nn_fct_h,   & 
     260            &                                                                      ' Vertical   order: ', nn_fct_v 
     261         CASE( np_FCT_zts )   ;   WRITE(numout,*) '         use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping' 
     262         CASE( np_MUS     )   ;   WRITE(numout,*) '         MUSCL    scheme is used' 
     263         CASE( np_UBS     )   ;   WRITE(numout,*) '         UBS      scheme is used' 
     264         CASE( np_QCK     )   ;   WRITE(numout,*) '         QUICKEST scheme is used' 
     265         END SELECT 
    266266      ENDIF 
    267267      ! 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90

    r5866 r5883  
    22   !!====================================================================== 
    33   !!                     ***  MODULE  traadv_cen  *** 
    4    !! Ocean  tracers:  horizontal & vertical advective trend (2nd/4th order centered) 
     4   !! Ocean  tracers:  advective trend (2nd/4th order centered) 
    55   !!====================================================================== 
    66   !! History :  3.7  ! 2014-05  (G. Madec)  original code 
     
    5252      !! ** Method  :   The advection is evaluated by a 2nd or 4th order scheme 
    5353      !!               using now fields (leap-frog scheme).  
    54       !! 
    5554      !!       kn_cen_h = 2  ==>> 2nd order centered scheme on the horizontal 
    5655      !!                = 4  ==>> 4th order    -        -       -      - 
    57       !! 
    5856      !!       kn_cen_v = 2  ==>> 2nd order centered scheme on the vertical 
    5957      !!                = 4  ==>> 4th order COMPACT  scheme     -      - 
    6058      !! 
    61       !! ** Action :  - update pta  with the now advective tracer trends 
    62       !!              - send trends to trdtra module for further diagnostcs 
     59      !! ** Action : - update pta  with the now advective tracer trends 
     60      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
     61      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
    6362      !!---------------------------------------------------------------------- 
    6463      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
     
    8988      ENDIF 
    9089      ! 
    91       !                          ! surface & bottom values  
    92       IF( .NOT.ln_linssh )   zwz(:,:, 1 ) = 0._wp             ! set to zero one for all 
    93                              zwz(:,:,jpk) = 0._wp             ! except at the surface in linear free surface 
     90      !                     
     91      zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers 
     92      zwz(:,:,jpk) = 0._wp 
    9493      ! 
    9594      DO jn = 1, kjpt            !==  loop over the tracers  ==! 
    9695         ! 
    97          SELECT CASE( kn_cen_h )          !--  Horizontal fluxes  --! 
    98          ! 
    99          CASE(  2  )                               ! 2nd order centered 
     96         SELECT CASE( kn_cen_h )       !--  Horizontal fluxes  --! 
     97         ! 
     98         CASE(  2  )                         !* 2nd order centered 
    10099            DO jk = 1, jpkm1 
    101100               DO jj = 1, jpjm1 
     
    107106            END DO 
    108107            ! 
    109          CASE(  4  )                               ! 4th order centered 
    110             ztu(:,:,jpk) = 0._wp                            ! Bottom value : flux set to zero 
     108         CASE(  4  )                         !* 4th order centered 
     109            ztu(:,:,jpk) = 0._wp                   ! Bottom value : flux set to zero 
    111110            ztv(:,:,jpk) = 0._wp 
    112             DO jk = 1, jpkm1                                 ! gradient 
    113                DO jj = 2, jpjm1                                   ! masked derivative 
     111            DO jk = 1, jpkm1                       ! masked gradient 
     112               DO jj = 2, jpjm1 
    114113                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    115114                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     
    120119            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn) 
    121120            ! 
    122             DO jk = 1, jpkm1                                 ! Horizontal advective fluxes 
     121            DO jk = 1, jpkm1                       ! Horizontal advective fluxes 
    123122               DO jj = 2, jpjm1 
    124123                  DO ji = 1, fs_jpim1   ! vector opt. 
     
    139138         END SELECT 
    140139         ! 
    141          !                             !==  Vertical fluxes  ==! 
    142          ! 
    143          SELECT CASE( kn_cen_v )             !* interior fluxes 
    144          ! 
    145          CASE(  2  )                               ! 2nd order centered 
     140         SELECT CASE( kn_cen_v )       !--  Vertical fluxes  --!   (interior) 
     141         ! 
     142         CASE(  2  )                         !* 2nd order centered 
    146143            DO jk = 2, jpk 
    147144               DO jj = 2, jpjm1 
     
    152149            END DO 
    153150            ! 
    154          CASE(  4  )                               ! 4th order centered 
    155             CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point 
     151         CASE(  4  )                         !* 4th order compact 
     152            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )      ! ztw = interpolated value of T at w-point 
    156153            DO jk = 2, jpkm1 
    157154               DO jj = 2, jpjm1 
     
    164161         END SELECT 
    165162         ! 
    166          IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask) 
     163         IF( ln_linssh ) THEN                !* top value   (linear free surf. only as zwz is multiplied by wmask) 
    167164            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean) 
    168165               DO jj = 1, jpj 
    169166                  DO ji = 1, jpi 
    170                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     167                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)  
    171168                  END DO 
    172169               END DO    
     
    182179                     &             - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )    & 
    183180                     &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )    & 
    184                      &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) / ( e1e2t(ji,jj) *  e3t_n(ji,jj,jk) ) 
     181                     &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    185182               END DO 
    186183            END DO 
    187184         END DO 
    188          !                                 ! trend diagnostics 
     185         !                             ! trend diagnostics 
    189186         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) THEN 
    190187            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
     
    192189            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 
    193190         END IF 
    194          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     191         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    195192         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    196193           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90

    r5866 r5883  
    5353      !!                  ***  ROUTINE tra_adv_fct  *** 
    5454      !!  
    55       !! **  Purpose :   Compute the now trend due to total advection of  
    56       !!       tracers and add it to the general trend of tracer equations 
     55      !! **  Purpose :   Compute the now trend due to total advection of tracers 
     56      !!               and add it to the general trend of tracer equations 
    5757      !! 
    5858      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction 
    5959      !!               (choice through the value of kn_fct) 
    60       !!               - 4th order compact scheme on the vertical  
     60      !!               - on the vertical the 4th order is a compact scheme  
    6161      !!               - corrected flux (monotonic correction)  
    6262      !! 
    63       !! ** Action : - update (pta) with the now advective tracer trends 
    64       !!             - send the trends for further diagnostics 
     63      !! ** Action : - update pta  with the now advective tracer trends 
     64      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
     65      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
    6566      !!---------------------------------------------------------------------- 
    6667      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
     
    101102      ENDIF 
    102103      ! 
    103       !                                         ! surface & bottom value : flux set to zero one for all 
    104       IF( .NOT.ln_linssh )   zwz(:,:, 1 ) = 0._wp                ! except at the surface in linear free surface case 
     104      !                          ! surface & bottom value : flux set to zero one for all 
     105      zwz(:,:, 1 ) = 0._wp             
    105106      zwx(:,:,jpk) = 0._wp   ;   zwy(:,:,jpk) = 0._wp    ;    zwz(:,:,jpk) = 0._wp 
    106107      ! 
    107108      zwi(:,:,:) = 0._wp         
    108       !                                                          ! =========== 
    109       DO jn = 1, kjpt                                            ! tracer loop 
    110          !                                                       ! =========== 
     109      ! 
     110      DO jn = 1, kjpt            !==  loop over the tracers  ==! 
    111111         ! 
    112112         !        !==  upstream advection with initial mass fluxes & intermediate update  ==! 
     
    126126         END DO 
    127127         !                    !* upstream tracer flux in the k direction *! 
    128          DO jk = 2, jpkm1         ! Interior value ( multiplied by wmask) 
     128         DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
    129129            DO jj = 1, jpj 
    130130               DO ji = 1, jpi 
     
    135135            END DO 
    136136         END DO 
    137          !                     
    138137         IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked) 
    139138            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface 
     
    155154                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    156155                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    157                      &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 
     156                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    158157                  ! update and guess with monotonic sheme 
    159158!!gm why tmask added in the two following lines ???    the mask is done in tranxt ! 
     
    174173         ENDIF 
    175174         ! 
    176          ! 
    177175         !        !==  anti-diffusive flux : high order minus low order  ==! 
    178176         ! 
    179          SELECT CASE( kn_fct_h )         !* horizontal anti-diffusive fluxes 
    180          ! 
    181          CASE(  2  )                         ! 2nd order centered 
     177         SELECT CASE( kn_fct_h )    !* horizontal anti-diffusive fluxes 
     178         ! 
     179         CASE(  2  )                   !- 2nd order centered 
    182180            DO jk = 1, jpkm1 
    183181               DO jj = 1, jpjm1 
     
    189187            END DO 
    190188            ! 
    191          CASE(  4  )                         ! 4th order centered 
    192             zltu(:,:,jpk) = 0._wp                            ! Bottom value : flux set to zero 
     189         CASE(  4  )                   !- 4th order centered 
     190            zltu(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero 
    193191            zltv(:,:,jpk) = 0._wp 
    194             DO jk = 1, jpkm1                                 ! Laplacian 
    195                DO jj = 1, jpjm1                                   ! First derivative (gradient) 
     192            DO jk = 1, jpkm1                 ! Laplacian 
     193               DO jj = 1, jpjm1                    ! 1st derivative (gradient) 
    196194                  DO ji = 1, fs_jpim1   ! vector opt. 
    197195                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     
    199197                  END DO 
    200198               END DO 
    201                DO jj = 2, jpjm1                                   !  
     199               DO jj = 2, jpjm1                    ! 2nd derivative * 1/ 6 
    202200                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    203201                     zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6 
     
    206204               END DO 
    207205            END DO 
    208             ! 
    209206            CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
    210207            ! 
    211             DO jk = 1, jpkm1                                 ! Horizontal advective fluxes 
     208            DO jk = 1, jpkm1                 ! Horizontal advective fluxes 
    212209               DO jj = 1, jpjm1 
    213210                  DO ji = 1, fs_jpim1   ! vector opt. 
     
    221218            END DO          
    222219            ! 
    223          CASE(  41 )                         ! 4th order centered       ==>>   !!gm coding attempt   need to be tested 
    224             ztu(:,:,jpk) = 0._wp                             ! Bottom value : flux set to zero 
     220         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested 
     221            ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero 
    225222            ztv(:,:,jpk) = 0._wp 
    226             DO jk = 1, jpkm1                                 ! gradient 
    227                DO jj = 1, jpjm1                                   ! First derivative (gradient) 
     223            DO jk = 1, jpkm1                 ! 1st derivative (gradient) 
     224               DO jj = 1, jpjm1 
    228225                  DO ji = 1, fs_jpim1   ! vector opt. 
    229226                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     
    234231            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn) 
    235232            ! 
    236             DO jk = 1, jpkm1                                 ! Horizontal advective fluxes 
     233            DO jk = 1, jpkm1                 ! Horizontal advective fluxes 
    237234               DO jj = 2, jpjm1 
    238235                  DO ji = 2, fs_jpim1   ! vector opt. 
     
    250247            ! 
    251248         END SELECT 
    252          !                                !* vertical anti-diffusive fluxes 
    253          SELECT CASE( kn_fct_v )                ! Interior values (w-masked) 
    254          ! 
    255          CASE(  2  )                                  ! 2nd order centered 
     249         !                       
     250         SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values) 
     251         ! 
     252         CASE(  2  )                   !- 2nd order centered 
    256253            DO jk = 2, jpkm1     
    257254               DO jj = 2, jpjm1 
    258255                  DO ji = fs_2, fs_jpim1 
    259                      zwz(ji,jj,jk) =  ( 0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   & 
    260                                        - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
    261                   END DO 
    262                END DO 
    263             END DO 
    264             ! 
    265          CASE(  4  )                                  ! 4th order COMPACT 
    266             !     
    267             CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! COMPACT interpolation of T at w-point 
    268             ! 
     256                     zwz(ji,jj,jk) =  (  pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   & 
     257                        &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk) 
     258                  END DO 
     259               END DO 
     260            END DO 
     261            ! 
     262         CASE(  4  )                   !- 4th order COMPACT 
     263            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )   ! zwt = COMPACT interpolation of T at w-point 
    269264            DO jk = 2, jpkm1 
    270265               DO jj = 2, jpjm1 
     
    276271            ! 
    277272         END SELECT 
    278          !                                      ! top ocean value: high order = upstream  ==>>  zwz=0 
    279          zwz(:,:, 1 ) = 0._wp                   ! only ocean surface as interior zwz values have been w-masked 
     273         IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0 
     274            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked 
     275         ENDIF 
    280276         ! 
    281277         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions 
    282278         CALL lbc_lnk( zwz, 'W',  1. ) 
    283  
     279         ! 
    284280         !        !==  monotonicity algorithm  ==! 
    285281         ! 
    286282         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 
    287  
    288  
     283         ! 
    289284         !        !==  final trend with corrected fluxes  ==! 
    290285         ! 
     
    300295         END DO 
    301296         ! 
    302          IF( l_trd ) THEN                 ! trend diagnostics (contribution of upstream fluxes) 
     297         IF( l_trd ) THEN     ! trend diagnostics (contribution of upstream fluxes) 
    303298            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed 
    304299            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
     
    311306            CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz ) 
    312307         END IF 
    313          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     308         !                    ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    314309         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    315            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:) 
    316            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:) 
     310           IF( jn == jp_tem )  htr_adv(:) = htr_adv(:) + ptr_sj( zwy(:,:,:) ) 
     311           IF( jn == jp_sal )  str_adv(:) = str_adv(:) + ptr_sj( zwy(:,:,:) ) 
    317312         ENDIF 
    318313         ! 
    319       END DO 
     314      END DO                     ! end of tracer loop 
    320315      ! 
    321316      CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 
     
    392387      zr_p2dt(:) = 1._wp / p2dt(:) 
    393388      ! 
     389      ! surface & Bottom value : flux set to zero for all tracers 
     390      zwz(:,:, 1 ) = 0._wp 
     391      zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp 
     392      zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp 
     393      ! 
    394394      !                                                          ! =========== 
    395395      DO jn = 1, kjpt                                            ! tracer loop 
    396396         !                                                       ! =========== 
    397          ! 1. Bottom value : flux set to zero 
    398          ! ---------------------------------- 
    399          zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp 
    400          zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp 
    401  
    402          ! 2. upstream advection with initial mass fluxes & intermediate update 
    403          ! -------------------------------------------------------------------- 
    404          ! upstream tracer flux in the i and j direction 
    405          DO jk = 1, jpkm1 
     397         ! 
     398         ! Upstream advection with initial mass fluxes & intermediate update 
     399         DO jk = 1, jpkm1        ! upstream tracer flux in the i and j direction 
    406400            DO jj = 1, jpjm1 
    407401               DO ji = 1, fs_jpim1   ! vector opt. 
     
    416410            END DO 
    417411         END DO 
    418  
    419          ! upstream tracer flux in the k direction 
    420          DO jk = 2, jpkm1         ! Interior value 
     412         !                       ! upstream tracer flux in the k direction 
     413         DO jk = 2, jpkm1              ! Interior value 
    421414            DO jj = 1, jpj 
    422415               DO ji = 1, jpi 
     
    427420            END DO 
    428421         END DO 
    429          !                       ! top value 
    430          IF( .NOT.ln_linssh ) THEN             ! variable volume: only k=1 as zwz is multiplied by wmask 
    431             zwz(:,:, 1 ) = 0._wp 
    432          ELSE                          ! linear free surface 
    433             IF( ln_isfcav ) THEN             ! ice-shelf cavities 
     422         IF( ln_linssh ) THEN          ! top value : linear free surface case only (as zwz is multiplied by wmask) 
     423            IF( ln_isfcav ) THEN             ! ice-shelf cavities: top value 
    434424               DO jj = 1, jpj 
    435425                  DO ji = 1, jpi 
     
    437427                  END DO 
    438428               END DO    
    439             ELSE                             ! standard case 
     429            ELSE                             ! no cavities, surface value 
    440430               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
    441431            ENDIF 
     
    446436            DO jj = 2, jpjm1 
    447437               DO ji = fs_2, fs_jpim1   ! vector opt. 
    448                   ! total intermediate advective trends 
     438                  !                             ! total intermediate advective trends 
    449439                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    450440                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    451                      &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 
    452                   ! update and guess with monotonic sheme 
     441                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     442                  !                             ! update and guess with monotonic sheme 
    453443                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra 
    454444                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 
     
    497487            END DO 
    498488         END DO 
    499        
     489         ! 
    500490         !                                !* vertical anti-diffusive flux 
    501491         zwz_sav(:,:,:)   = zwz(:,:,:) 
    502492         ztrs   (:,:,:,1) = ptb(:,:,:,jn) 
    503493         zwzts  (:,:,:)   = 0._wp 
    504          IF( .NOT.ln_linssh )   zwz(:,:, 1 ) = 0._wp    ! surface value set to zero in vvl case 
    505494         ! 
    506495         DO jl = 1, kn_fct_zts                  ! Start of sub timestepping loop 
     
    535524                     END DO 
    536525                  END DO    
    537                ELSE                                      ! standard case 
     526               ELSE                                      ! no ocean cavities 
    538527                  zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
    539528               ENDIF 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90

    r5866 r5883  
    6262      !!              ld_msc_ups=T :  
    6363      !! 
    64       !! ** Action  : - update (ta,sa) with the now advective tracer trends 
    65       !!              - send trends to trdtra module for further diagnostcs 
     64      !! ** Action : - update pta  with the now advective tracer trends 
     65      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
     66      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
    6667      !! 
    6768      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation 
     
    116117      ENDIF  
    117118      !       
    118       !                                                     ! =========== 
    119       DO jn = 1, kjpt                                       ! tracer loop 
    120          !                                                  ! =========== 
    121          ! I. Horizontal advective fluxes 
    122          ! ------------------------------ 
    123          ! first guess of the slopes 
    124          zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0        ! bottom values 
    125          ! interior values 
    126          DO jk = 1, jpkm1 
     119      DO jn = 1, kjpt            !==  loop over the tracers  ==! 
     120         ! 
     121         !                          !* Horizontal advective fluxes 
     122         ! 
     123         !                                !-- first guess of the slopes 
     124         zwx(:,:,jpk) = 0.e0                    ! bottom values 
     125         zwy(:,:,jpk) = 0._wp   
     126         DO jk = 1, jpkm1                       ! interior values 
    127127            DO jj = 1, jpjm1       
    128128               DO ji = 1, fs_jpim1   ! vector opt. 
     
    132132           END DO 
    133133         END DO 
    134          ! 
    135          CALL lbc_lnk( zwx, 'U', -1. )                        ! lateral boundary conditions on zwx, zwy   (changed sign) 
     134         CALL lbc_lnk( zwx, 'U', -1. )          ! lateral boundary conditions   (changed sign) 
    136135         CALL lbc_lnk( zwy, 'V', -1. ) 
    137          !                                             !-- Slopes of tracer 
    138          zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0    ! bottom values 
    139          DO jk = 1, jpkm1                                     ! interior values 
     136         !                                !-- Slopes of tracer 
     137         zslpx(:,:,jpk) = 0._wp                  ! bottom values 
     138         zslpy(:,:,jpk) = 0._wp 
     139         DO jk = 1, jpkm1                        ! interior values 
    140140            DO jj = 2, jpj 
    141141               DO ji = fs_2, jpi   ! vector opt. 
     
    148148         END DO 
    149149         ! 
    150          DO jk = 1, jpkm1                                     ! Slopes limitation 
     150         DO jk = 1, jpkm1                 !-- Slopes limitation 
    151151            DO jj = 2, jpj 
    152152               DO ji = fs_2, jpi   ! vector opt. 
     
    161161         END DO 
    162162         ! 
    163          !                                             !-- MUSCL horizontal advective fluxes 
    164          DO jk = 1, jpkm1                                     ! interior values 
     163         DO jk = 1, jpkm1                 !-- MUSCL horizontal advective fluxes 
    165164            zdt  = p2dt(jk) 
    166165            DO jj = 2, jpjm1 
     
    185184         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign) 
    186185         ! 
    187          DO jk = 1, jpkm1         ! Tracer flux divergence at t-point added to the general trend 
     186         DO jk = 1, jpkm1                 !-- Tracer advective trend 
    188187            DO jj = 2, jpjm1       
    189188               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    194193           END DO 
    195194         END DO         
    196          !                                 ! trend diagnostics (contribution of upstream fluxes) 
     195         !                                ! trend diagnostics 
    197196         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   & 
    198197            &( cdtype == 'TRC' .AND. l_trdtrc )      )  THEN 
     
    200199            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 
    201200         END IF 
    202          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     201         !                                ! "Poleward" heat and salt transports 
    203202         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    204203            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
    205204            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    206205         ENDIF 
    207  
    208          ! II. Vertical advective fluxes 
    209          ! ----------------------------- 
     206         ! 
     207         !                          !* Vertical advective fluxes 
     208         ! 
    210209         !                                !-- first guess of the slopes 
    211210         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions 
    212          zwx(:,:,jpk) = 0._wp                   ! surface & bottom boundary conditions 
    213          DO jk = 2, jpkm1                                     ! interior values 
     211         zwx(:,:,jpk) = 0._wp 
     212         DO jk = 2, jpkm1                       ! interior values 
    214213            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) 
    215214         END DO 
    216  
    217215         !                                !-- Slopes of tracer 
    218216         zslpx(:,:,1) = 0._wp                   ! surface values 
     
    220218            DO jj = 1, jpj 
    221219               DO ji = 1, jpi 
    222                   zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   & 
    223                      &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 
    224                END DO 
    225             END DO 
    226          END DO 
    227          !                                !-- Slopes limitation 
    228          DO jk = 2, jpkm1                       ! interior values 
    229             DO jj = 1, jpj 
     220                  zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  & 
     221                     &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  ) 
     222               END DO 
     223            END DO 
     224         END DO 
     225         DO jk = 2, jpkm1                 !-- Slopes limitation 
     226            DO jj = 1, jpj                      ! interior values 
    230227               DO ji = 1, jpi 
    231228                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
     
    235232            END DO 
    236233         END DO 
    237          !                                !-- vertical advective flux 
    238          DO jk = 1, jpkm1                       ! interior values 
     234         DO jk = 1, jpk-2                 !-- vertical advective flux 
    239235            zdt  = p2dt(jk) 
    240236            DO jj = 2, jpjm1       
     
    242238                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 
    243239                  zalpha = 0.5 + z0w 
    244                   zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt  / ( e1e2t(ji,jj) * e3w_n(ji,jj,jk+1) ) 
     240                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1) 
    245241                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
    246242                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
     
    249245            END DO 
    250246         END DO 
    251          !                                      ! top values  (bottom already set to zero) 
    252          IF( ln_linssh ) THEN                         !* linear free surface  
    253             IF( ln_isfcav ) THEN                            ! ice-shelf cavities (top of the ocean) 
     247         IF( ln_linssh ) THEN                   ! top values, linear free surface only 
     248            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean) 
    254249               DO jj = 1, jpj 
    255250                  DO ji = 1, jpi 
     
    257252                  END DO 
    258253               END DO    
    259             ELSE                                             ! no cavities: only at the ocean surface 
     254            ELSE                                      ! no cavities: only at the ocean surface 
    260255               zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
    261256            ENDIF 
    262          ELSE                                       !* non-linear free surface 
    263             zwx(:,:, 1 ) = 0._wp                            ! k=1 only as zwx has been multiplied by wmask 
    264257         ENDIF 
    265258         ! 
    266          DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend 
     259         DO jk = 1, jpkm1                 !-- vertical advective trend 
    267260            DO jj = 2, jpjm1       
    268261               DO ji = fs_2, fs_jpim1   ! vector opt. 
    269                   pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 
    270                END DO 
    271             END DO 
    272          END DO 
    273          !                                 ! Save the vertical advective trends for diagnostic 
     262                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     263               END DO 
     264            END DO 
     265         END DO 
     266         !                                ! send trends for diagnostic 
    274267         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     & 
    275268            &( cdtype == 'TRC' .AND. l_trdtrc )      )   & 
    276269            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 
    277270         ! 
    278       END DO 
     271      END DO                     ! end of tracer loop 
    279272      ! 
    280273      CALL wrk_dealloc( jpi,jpj,jpk,   zslpx, zslpy, zwx, zwy ) 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90

    r5866 r5883  
    7878      !!            prevent the appearance of spurious numerical oscillations 
    7979      !! 
    80       !! ** Action : - update (pta) with the now advective tracer trends 
    81       !!             - save the trends  
     80      !! ** Action : - update pta  with the now advective tracer trends 
     81      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
     82      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
    8283      !! 
    8384      !! ** Reference : Leonard (1979, 1991) 
     
    105106      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
    106107      ! 
    107       ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 
     108      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 
    108109      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt )  
    109110      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt )  
    110111 
    111       ! II. The vertical fluxes are computed with the 2nd order centered scheme 
     112      !        ! vertical fluxes are computed with the 2nd order centered scheme 
    112113      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt ) 
    113114      ! 
     
    224225            END DO 
    225226         END DO 
    226          !                                 ! trend diagnostics (contribution of upstream fluxes) 
     227         !                                 ! trend diagnostics 
    227228         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
    228229         ! 
     
    348349            END DO 
    349350         END DO 
    350          !                                 ! trend diagnostics (contribution of upstream fluxes) 
     351         !                                 ! trend diagnostics 
    351352         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
    352353         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     
    381382      CALL wrk_alloc( jpi,jpj,jpk,   zwz ) 
    382383      ! 
    383       !                          ! surface & bottom values  
    384       IF( .NOT.ln_linssh )   zwz(:,:, 1 ) = 0._wp             ! set to zero one for all 
    385                              zwz(:,:,jpk) = 0._wp             ! except at the surface in linear free surface 
     384      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers 
     385      zwz(:,:,jpk) = 0._wp 
    386386      ! 
    387387      !                                                          ! =========== 
     
    403403                  END DO 
    404404               END DO    
    405             ELSE                                   ! no ice-shelf cavities (only ocean surface) 
     405            ELSE                                   ! no ocean cavities (only ocean surface) 
    406406               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 
    407407            ENDIF 
     
    416416            END DO 
    417417         END DO 
    418          !                                 ! Save the vertical advective trends for diagnostic 
     418         !                                 ! Send trends for diagnostic 
    419419         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 
    420420         ! 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

    r5866 r5883  
    7171      !!                On the vertical, the advection is evaluated using a FCT scheme, 
    7272      !!      as the UBS have been found to be too diffusive.  
    73 !!gm  !!                kn_ubs_v argument (not coded for the moment) 
    74       !!      controles whether the FCT is based on a 2nd order centrered scheme (kn_ubs_v=2)  
    75       !!      or on a 4th order compact scheme (kn_ubs_v=4). 
     73      !!                kn_ubs_v argument controles whether the FCT is based on  
     74      !!      a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact  
     75      !!      scheme (kn_ubs_v=4). 
    7676      !! 
    77       !! ** Action : - update (pta) with the now advective tracer trends 
     77      !! ** Action : - update pta  with the now advective tracer trends 
     78      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
     79      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
    7880      !! 
    7981      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.  
     
    110112      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    111113      ! 
    112       zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp     ! Bottom value : set to zero one for all 
     114      ztw (:,:, 1 ) = 0._wp      ! surface & bottom value : set to zero for all tracers 
     115      zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp 
    113116      ztw (:,:,jpk) = 0._wp   ;   zti (:,:,jpk) = 0._wp 
    114       IF( .NOT.ln_linssh )   ztw(:,:, 1 ) = 0._wp                   ! surface value: set to zero only in vvl case 
    115117      ! 
    116118      !                                                          ! =========== 
     
    264266                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk)                          & 
    265267                        &           + ptn(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   & 
    266                         &                              / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 
     268                        &                              * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    267269                  END DO 
    268270               END DO 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

    r5845 r5883  
    1212 
    1313   !!---------------------------------------------------------------------- 
    14    !!   tra_bbc      : update the tracer trend at ocean bottom  
    15    !!   tra_bbc_init : initialization of geothermal heat flux trend 
     14   !!   tra_bbc       : update the tracer trend at ocean bottom  
     15   !!   tra_bbc_init  : initialization of geothermal heat flux trend 
    1616   !!---------------------------------------------------------------------- 
    17    USE oce             ! ocean variables 
    18    USE dom_oce         ! domain: ocean 
    19    USE phycst          ! physical constants 
    20    USE trd_oce         ! trends: ocean variables 
    21    USE trdtra          ! trends manager: tracers  
    22    USE in_out_manager  ! I/O manager 
    23    USE iom             ! I/O manager 
    24    USE fldread         ! read input fields 
    25    USE lbclnk            ! ocean lateral boundary conditions (or mpp link) 
    26    USE lib_mpp           ! distributed memory computing library 
    27    USE prtctl          ! Print control 
    28    USE wrk_nemo        ! Memory Allocation 
    29    USE timing          ! Timing 
     17   USE oce            ! ocean variables 
     18   USE dom_oce        ! domain: ocean 
     19   USE phycst         ! physical constants 
     20   USE trd_oce        ! trends: ocean variables 
     21   USE trdtra         ! trends manager: tracers  
     22   ! 
     23   USE in_out_manager ! I/O manager 
     24   USE iom            ! xIOS  
     25   USE fldread        ! read input fields 
     26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     27   USE lib_mpp        ! distributed memory computing library 
     28   USE prtctl         ! Print control 
     29   USE wrk_nemo       ! Memory Allocation 
     30   USE timing         ! Timing 
    3031 
    3132   IMPLICIT NONE 
     
    4243   REAL(wp), PUBLIC , ALLOCATABLE, DIMENSION(:,:) ::   qgh_trd0   ! geothermal heating trend 
    4344 
    44    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh              ! structure of input qgh (file informations, fields read) 
     45   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh   ! structure of input qgh (file informations, fields read) 
    4546  
    4647   !!---------------------------------------------------------------------- 
     
    6768      !!       Where Qsf is the geothermal heat flux. 
    6869      !! 
    69       !! ** Action  : - update the temperature trends (ta) with the trend of 
    70       !!                the ocean bottom boundary condition 
     70      !! ** Action  : - update the temperature trends with geothermal heating trend 
     71      !!              - send the trend for further diagnostics (ln_trdtra=T) 
    7172      !! 
    7273      !! References : Stein, C. A., and S. Stein, 1992, Nature, 359, 123-129. 
     
    7475      !!---------------------------------------------------------------------- 
    7576      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    76       !! 
    77       INTEGER  ::   ji, jj, ik    ! dummy loop indices 
    78       REAL(wp) ::   zqgh_trd      ! geothermal heat flux trend 
     77      ! 
     78      INTEGER  ::   ji, jj    ! dummy loop indices 
    7979      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt 
    8080      !!---------------------------------------------------------------------- 
     
    8282      IF( nn_timing == 1 )  CALL timing_start('tra_bbc') 
    8383      ! 
    84       IF( l_trdtra )   THEN         ! Save ta and sa trends 
    85          CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 
     84      IF( l_trdtra )   THEN         ! Save the input temperature trend 
     85         CALL wrk_alloc( jpi,jpj,jpk,  ztrdt ) 
    8686         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    8787      ENDIF 
    88       ! 
    89       !                             !  Add the geothermal heat flux trend on temperature 
     88      !                             !  Add the geothermal trend on temperature 
    9089      DO jj = 2, jpjm1 
    9190         DO ji = 2, jpim1 
    92             ik = mbkt(ji,jj) 
    93             zqgh_trd = qgh_trd0(ji,jj) / e3t_n(ji,jj,ik) 
    94             tsa(ji,jj,ik,jp_tem) = tsa(ji,jj,ik,jp_tem) + zqgh_trd 
     91            tsa(ji,jj,mbkt(ji,jj),jp_tem) = tsa(ji,jj,mbkt(ji,jj),jp_tem) + qgh_trd0(ji,jj) / e3t_n(ji,jj,mbkt(ji,jj)) 
    9592         END DO 
    9693      END DO 
     
    9895      CALL lbc_lnk( tsa(:,:,:,jp_tem) , 'T', 1. ) 
    9996      ! 
    100       IF( l_trdtra ) THEN        ! Save the geothermal heat flux trend for diagnostics 
     97      IF( l_trdtra ) THEN        ! Send the trend for diagnostics 
    10198         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    10299         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt ) 
    103          CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 
     100         CALL wrk_dealloc( jpi,jpj,jpk,  ztrdt ) 
    104101      ENDIF 
    105102      ! 
     
    126123      !! ** Action  : - read/fix the geothermal heat qgh_trd0 
    127124      !!---------------------------------------------------------------------- 
    128       USE iom 
    129       !! 
    130125      INTEGER  ::   ji, jj              ! dummy loop indices 
    131126      INTEGER  ::   inum                ! temporary logical unit 
     
    138133      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir  
    139134      !!---------------------------------------------------------------------- 
    140  
     135      ! 
    141136      REWIND( numnam_ref )              ! Namelist nambbc in reference namelist : Bottom momentum boundary condition 
    142137      READ  ( numnam_ref, nambbc, IOSTAT = ios, ERR = 901) 
    143138901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbc in reference namelist', lwp ) 
    144  
     139      ! 
    145140      REWIND( numnam_cfg )              ! Namelist nambbc in configuration namelist : Bottom momentum boundary condition 
    146141      READ  ( numnam_cfg, nambbc, IOSTAT = ios, ERR = 902 ) 
    147142902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbc in configuration namelist', lwp ) 
    148143      IF(lwm) WRITE ( numond, nambbc ) 
    149  
     144      ! 
    150145      IF(lwp) THEN                     ! Control print 
    151146         WRITE(numout,*) 
     
    158153         WRITE(numout,*) 
    159154      ENDIF 
    160  
     155      ! 
    161156      IF( ln_trabbc ) THEN             !==  geothermal heating  ==! 
    162157         ! 
     
    189184            WRITE(ctmp1,*) '     bad flag value for nn_geoflx = ', nn_geoflx 
    190185            CALL ctl_stop( ctmp1 ) 
    191             ! 
    192186         END SELECT 
    193187         ! 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r5845 r5883  
    111111      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl') 
    112112      ! 
    113       IF( l_trdtra )   THEN                         !* Save ta and sa trends 
     113      IF( l_trdtra )   THEN                         !* Save the input trends 
    114114         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    115115         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
     
    131131         ! 
    132132      END IF 
    133  
     133      ! 
    134134      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl 
    135135         ! 
     
    145145      END IF 
    146146 
    147       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
     147      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics 
    148148         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    149149         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
     
    301301            ! 
    302302         END DO 
    303          !                                                  ! =========== 
    304       END DO                                                ! end tracer 
    305       !                                                     ! =========== 
    306       ! 
     303         !                                                       ! =========== 
     304      END DO                                                     ! end tracer 
     305      !                                                          ! =========== 
    307306      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv') 
    308307      ! 
     
    339338      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index 
    340339      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    341       !! 
     340      ! 
    342341      INTEGER  ::   ji, jj                    ! dummy loop indices 
    343342      INTEGER  ::   ik                        ! local integers 
     
    400399         ! 
    401400      ENDIF 
    402  
     401      ! 
    403402      !                                   !-------------------! 
    404403      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   ! 
     
    499498      INTEGER ::   ios                  !   -      - 
    500499      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 
    501       !! 
     500      ! 
    502501      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl 
    503502      !!---------------------------------------------------------------------- 
     
    505504      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init') 
    506505      ! 
    507       CALL wrk_alloc( jpi, jpj, zmbk ) 
    508       ! 
    509  
    510506      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme 
    511507      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901) 
     
    544540      END DO 
    545541      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
     542      CALL wrk_alloc( jpi, jpj, zmbk ) 
    546543      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    547544      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
     545      CALL wrk_dealloc( jpi, jpj, zmbk ) 
    548546 
    549547                                        !* sign of grad(H) at u- and v-points 
     
    592590      ENDIF 
    593591      ! 
    594       CALL wrk_dealloc( jpi, jpj, zmbk ) 
    595       ! 
    596592      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init') 
    597593      ! 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r5845 r5883  
    3131   USE dtatsd         ! data: temperature & salinity 
    3232   USE zdfmxl         ! vertical physics: mixed layer depth 
     33   ! 
    3334   USE in_out_manager ! I/O manager 
    3435   USE lib_mpp        ! MPP library 
     
    4142   PRIVATE 
    4243 
    43    PUBLIC   tra_dmp      ! routine called by step.F90 
    44    PUBLIC   tra_dmp_init ! routine called by nemogcm.F90 
     44   PUBLIC   tra_dmp        ! called by step.F90 
     45   PUBLIC   tra_dmp_init   ! called by nemogcm.F90 
    4546 
    4647   !                                           !!* Namelist namtra_dmp : T & S newtonian damping * 
     
    8889      !!      below the well mixed layer (nlmdmp=2) 
    8990      !! 
    90       !! ** Action  : - (ta,sa)  tracer trends updated with the damping trend 
     91      !! ** Action  : - tsa: tracer trends updated with the damping trend 
    9192      !!---------------------------------------------------------------------- 
    9293      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     
    176177      !!---------------------------------------------------------------------- 
    177178      INTEGER ::   ios, imask   ! local integers  
    178       !! 
     179      ! 
    179180      NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto 
    180181      !!---------------------------------------------------------------------- 
     
    228229   END SUBROUTINE tra_dmp_init 
    229230 
     231   !!====================================================================== 
    230232END MODULE tradmp 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r5845 r5883  
    1212 
    1313   !!---------------------------------------------------------------------- 
    14    !!   tra_ldf      : update the tracer trend with the lateral diffusion 
    15    !!   tra_ldf_init : initialization, namelist read, and parameters control 
    16    !!---------------------------------------------------------------------- 
    17    USE oce           ! ocean dynamics and tracers 
    18    USE dom_oce       ! ocean space and time domain 
    19    USE phycst        ! physical constants 
    20    USE ldftra        ! lateral diffusion: eddy diffusivity & EIV coeff. 
    21    USE ldfslp        ! lateral diffusion: iso-neutral slope 
    22    USE traldf_lap    ! lateral diffusion: laplacian iso-level            operator  (tra_ldf_lap   routine) 
    23    USE traldf_iso    ! lateral diffusion: laplacian iso-neutral standard operator  (tra_ldf_iso   routine) 
    24    USE traldf_triad  ! lateral diffusion: laplacian iso-neutral triad    operator  (tra_ldf_triad routine) 
    25    USE traldf_blp    ! lateral diffusion (iso-level lap/blp)                       (tra_ldf_lap   routine) 
    26    USE trd_oce       ! trends: ocean variables 
    27    USE trdtra        ! ocean active tracers trends 
     14   !!   tra_ldf       : update the tracer trend with the lateral diffusion trend 
     15   !!   tra_ldf_init  : initialization, namelist read, and parameters control 
     16   !!---------------------------------------------------------------------- 
     17   USE oce            ! ocean dynamics and tracers 
     18   USE dom_oce        ! ocean space and time domain 
     19   USE phycst         ! physical constants 
     20   USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff. 
     21   USE ldfslp         ! lateral diffusion: iso-neutral slope 
     22   USE traldf_lap_blp ! lateral diffusion: laplacian iso-level            operator  (tra_ldf_lap/_blp   routines) 
     23   USE traldf_iso     ! lateral diffusion: laplacian iso-neutral standard operator  (tra_ldf_iso        routine ) 
     24   USE traldf_triad   ! lateral diffusion: laplacian iso-neutral triad    operator  (tra_ldf_triad      routine ) 
     25   USE trd_oce        ! trends: ocean variables 
     26   USE trdtra         ! ocean active tracers trends 
    2827   ! 
    2928   USE prtctl         ! Print control 
     
    7170      ! 
    7271      SELECT CASE ( nldf )                     !* compute lateral mixing trend and add it to the general trend 
    73       ! 
    7472      CASE ( np_lap   )                                  ! laplacian: iso-level operator 
    7573         CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   ) 
     
    8179         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf ) 
    8280      END SELECT 
    83  
     81      ! 
    8482      IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics 
    8583         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
     
    113111         WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator' 
    114112         WRITE(numout,*) '~~~~~~~~~~~' 
    115          WRITE(numout,*) '   Namelist namtra_ldf already read in ldftra module' 
    116          WRITE(numout,*) '   see ldf_tra_init report for lateral mixing parameters' 
     113         WRITE(numout,*) '   Namelist namtra_ldf: already read in ldftra module' 
     114         WRITE(numout,*) '      see ldf_tra_init report for lateral mixing parameters' 
    117115         WRITE(numout,*) 
    118116      ENDIF 
     
    177175      ENDIF 
    178176      ! 
    179       IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' ) 
     177      IF( ierr == 1 )   CALL ctl_stop( 'iso-level in z-partial step, not allowed' ) 
    180178      IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) )                                    & 
    181            &            CALL ctl_stop( '          eddy induced velocity on tracers requires isopycnal',    & 
    182            &                                                                    ' laplacian diffusion' ) 
     179           &            CALL ctl_stop( 'eddy induced velocity on tracers requires iso-neutral laplacian diffusion' ) 
     180           ! 
    183181      IF(  nldf == np_lap_i .OR. nldf == np_lap_it .OR. & 
    184182         & nldf == np_blp_i .OR. nldf == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required  
     
    186184      IF(lwp) THEN 
    187185         WRITE(numout,*) 
    188          IF( nldf == np_no_ldf )   WRITE(numout,*) '          NO lateral diffusion' 
    189          IF( nldf == np_lap    )   WRITE(numout,*) '          laplacian iso-level operator' 
    190          IF( nldf == np_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)' 
    191          IF( nldf == np_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)' 
    192          IF( nldf == np_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator' 
    193          IF( nldf == np_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)' 
    194          IF( nldf == np_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)' 
     186         SELECT CASE( nldf ) 
     187         CASE( np_no_ldf )   ;   WRITE(numout,*) '   NO lateral diffusion' 
     188         CASE( np_lap    )   ;   WRITE(numout,*) '   laplacian iso-level operator' 
     189         CASE( np_lap_i  )   ;   WRITE(numout,*) '   Rotated laplacian operator (standard)' 
     190         CASE( np_lap_it )   ;   WRITE(numout,*) '   Rotated laplacian operator (triad)' 
     191         CASE( np_blp    )   ;   WRITE(numout,*) '   bilaplacian iso-level operator' 
     192         CASE( np_blp_i  )   ;   WRITE(numout,*) '   Rotated bilaplacian operator (standard)' 
     193         CASE( np_blp_it )   ;   WRITE(numout,*) '   Rotated bilaplacian operator (triad)' 
     194         END SELECT 
    195195      ENDIF 
    196196      ! 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r5845 r5883  
    1414 
    1515   !!---------------------------------------------------------------------- 
    16    !!   tra_ldf_iso  : update the tracer trend with the horizontal component of a iso-neutral laplacian operator 
    17    !!                  and with the vertical part of the isopycnal or geopotential s-coord. operator  
     16   !!   tra_ldf_iso   : update the tracer trend with the horizontal component of a iso-neutral laplacian operator 
     17   !!                   and with the vertical part of the isopycnal or geopotential s-coord. operator  
    1818   !!---------------------------------------------------------------------- 
    19    USE oce             ! ocean dynamics and active tracers 
    20    USE dom_oce         ! ocean space and time domain 
    21    USE trc_oce         ! share passive tracers/Ocean variables 
    22    USE zdf_oce         ! ocean vertical physics 
    23    USE ldftra          ! lateral diffusion: tracer eddy coefficients 
    24    USE ldfslp          ! iso-neutral slopes 
    25    USE diaptr          ! poleward transport diagnostics 
     19   USE oce            ! ocean dynamics and active tracers 
     20   USE dom_oce        ! ocean space and time domain 
     21   USE trc_oce        ! share passive tracers/Ocean variables 
     22   USE zdf_oce        ! ocean vertical physics 
     23   USE ldftra         ! lateral diffusion: tracer eddy coefficients 
     24   USE ldfslp         ! iso-neutral slopes 
     25   USE diaptr         ! poleward transport diagnostics 
    2626   ! 
    27    USE in_out_manager  ! I/O manager 
    28    USE iom             ! I/O library 
    29    USE phycst          ! physical constants 
    30    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    31    USE wrk_nemo        ! Memory Allocation 
    32    USE timing          ! Timing 
     27   USE in_out_manager ! I/O manager 
     28   USE iom            ! I/O library 
     29   USE phycst         ! physical constants 
     30   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     31   USE wrk_nemo       ! Memory Allocation 
     32   USE timing         ! Timing 
    3333 
    3434   IMPLICIT NONE 
     
    126126         ah_wslp2(:,:,:) = 0._wp 
    127127      ENDIF 
    128       ! 
    129128      !                                               ! set time step size (Euler/Leapfrog) 
    130129      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdttra(1)      ! at nit000   (Euler) 
     
    136135      ELSE                    ;   zsign = -1._wp 
    137136      ENDIF 
    138           
    139137          
    140138      !!---------------------------------------------------------------------- 
     
    241239            ENDIF 
    242240         ENDIF 
    243  
     241         ! 
    244242         !!---------------------------------------------------------------------- 
    245243         !!   II - horizontal trend  (full) 
     
    265263!         END IF 
    266264!!gm  
    267  
    268265            DO jj = 1 , jpjm1            !==  Horizontal fluxes 
    269266               DO ji = 1, fs_jpim1   ! vector opt. 
     
    298295         END DO                                        !   End of slab   
    299296 
    300  
    301297         !!---------------------------------------------------------------------- 
    302298         !!   III - vertical trend (full) 
    303299         !!---------------------------------------------------------------------- 
    304           
     300         ! 
    305301         ztfw(1,:,:) = 0._wp     ;     ztfw(jpi,:,:) = 0._wp 
    306           
     302         ! 
    307303         ! Vertical fluxes 
    308304         ! --------------- 
    309           
    310          ! Surface and bottom vertical fluxes set to zero 
     305         !                          ! Surface and bottom vertical fluxes set to zero 
    311306         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp 
    312307          
    313          ! interior (2=<jk=<jpk-1) 
    314          DO jk = 2, jpkm1 
     308         DO jk = 2, jpkm1           ! interior (2=<jk=<jpk-1) 
    315309            DO jj = 2, jpjm1 
    316310               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    336330            END DO 
    337331         END DO 
    338          ! 
    339332         !                                !==  add the vertical 33 flux  ==! 
    340333         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90

    </
    r5861 r5883  
    1 MODULE traldf_lap 
     1MODULE traldf_lap_blp 
    22   !!============================================================================== 
    3    !!                       ***  MODULE  traldf_lap  *** 
     3   !!                       ***  MODULE  traldf_lap_blp  *** 
    44   !! Ocean tracers:  lateral diffusivity trend  (laplacian and bilaplacian) 
    55   !!============================================================================== 
    6    !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code 
    7    !!                 ! 1991-11  (G. Madec) 
    8    !!                 ! 1995-11  (G. Madec)  suppress volumetric scale factors 
    9    !!                 ! 1996-01  (G. Madec)  statement function for e3 
    10    !!            NEMO ! 2002-06  (G. Madec)  F90: Free form and module 
    11    !!            1.0  ! 2004-08  (C. Talandier) New trends organization 
    12    !!                 ! 2005-11  (G. Madec)  add zps case 
    13    !!            3.0  ! 2010-06  (C. Ethe, G. Madec) Merge TRA-TRC 
    14    !!            3.7  ! 2014-01  (G. Madec, S. Masson) re-entrant laplacian  
    15    !!---------------------------------------------------------------------- 
    16  
    17    !!---------------------------------------------------------------------- 
    18    !!   tra_ldf_lap : update the tracer trend with the lateral diffusion : iso-level laplacian operator 
    19    !!   tra_ldf_blp : update the tracer trend with the lateral diffusion : iso-level bilaplacian operator 
    20    !!---------------------------------------------------------------------- 
    21    USE oce             ! ocean dynamics and active tracers 
    22    USE dom_oce         ! ocean space and time domain 
    23    USE ldftra          ! lateral physics: eddy diffusivity 
    24    USE diaptr          ! poleward transport diagnostics 
    25    USE trc_oce         ! share passive tracers/Ocean variables 
    26    USE zpshde          ! partial step: hor. derivative     (zps_hde routine) 
     6   !! History :  3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian  
     7   !!---------------------------------------------------------------------- 
     8 
     9   !!---------------------------------------------------------------------- 
     10   !!   tra_ldf_lap   : tracer trend update with iso-level laplacian diffusive operator 
     11   !!   tra_ldf_blp   : tracer trend update with iso-level or iso-neutral bilaplacian operator 
     12   !!---------------------------------------------------------------------- 
     13   USE oce            ! ocean dynamics and active tracers 
     14   USE dom_oce        ! ocean space and time domain 
     15   USE ldftra         ! lateral physics: eddy diffusivity 
     16   USE traldf_iso     ! iso-neutral lateral diffusion (standard operator)     (tra_ldf_iso   routine)