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

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

Changeset 9017


Ignore:
Timestamp:
2017-12-13T15:37:41+01:00 (4 years ago)
Author:
davestorkey
Message:

UKMO/dev_r8600_closea_rewrite branch : commit code

Location:
branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM
Files:
1 added
4 edited
1 moved

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

    r8995 r9017  
    1 MODULE usrdef_closea 
     1MODULE closea 
    22   !!====================================================================== 
    33   !!                   ***  MODULE  usrdef_closea  *** 
    4    !! 
    5    !!                      ===  ORCA configuration  === 
    6    !!                         (2, 1 and 1/4 degrees) 
    74   !! 
    85   !! User define : specific treatments associated with closed seas 
     
    1310   !!             3.4  !  2014-12  (P.G. Fogli) sbc_clo bug fix & mpp reproducibility 
    1411   !!             4.0  !  2016-06  (G. Madec)  move to usrdef_closea, remove clo_ups 
     12   !!             4.0  !  2017-12  (D. Storkey) new formulation based on masks read from file 
    1513   !!---------------------------------------------------------------------- 
    1614 
     
    2523   USE phycst          ! physical constants 
    2624   USE sbc_oce         ! ocean surface boundary conditions 
     25   USE iom             ! I/O routines 
    2726   ! 
    2827   USE in_out_manager  ! I/O manager 
    29    USE lib_fortran,    ONLY: glob_sum, DDPDD 
     28   USE lib_fortran,    ONLY: glob_sum 
    3029   USE lbclnk          ! lateral boundary condition - MPP exchanges 
    3130   USE lib_mpp         ! MPP library 
    3231   USE timing 
     32   USE wrk_nemo        ! Memory allocation 
    3333 
    3434   IMPLICIT NONE 
     
    3636 
    3737   PUBLIC dom_clo      ! called by domain module 
    38    PUBLIC sbc_clo      ! called by step module 
     38   PUBLIC sbc_clo      ! called by sbcmod module 
    3939   PUBLIC clo_rnf      ! called by sbcrnf module 
    40    PUBLIC clo_bat      ! called in domzgr module 
    41  
    42    INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 4      !: number of closed sea 
    43    INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea 
    44    INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j) 
    45    INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi2, ncsj2     !: north-east closed sea limits (i,j) 
    46    INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsnr            !: number of point where run-off pours 
    47    INTEGER, PUBLIC, DIMENSION(jpncs,4) ::   ncsir, ncsjr     !: Location of runoff 
    48  
    49    REAL(wp), DIMENSION (jpncs+1)       ::   surf             ! closed sea surface 
     40   PUBLIC clo_bat      ! called in domain module 
     41 
     42   INTEGER, PUBLIC :: jncs !: number of closed seas (inferred from closea_mask field) 
     43   INTEGER, PUBLIC :: jncsr !: number of closed seas rnf mappings (inferred from closea_mask_rnf field) 
     44   INTEGER, PUBLIC :: jncse !: number of closed seas empmr mappings (inferred from closea_mask_empmr field) 
     45    
     46   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask       !: mask of integers defining closed seas 
     47   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask_rnf   !: mask of integers defining closed seas rnf mappings 
     48   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask_empmr !: mask of integers defining closed seas empmr mappings 
     49   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surf         !: closed sea surface areas  
     50                                                                  !: (and residual global surface area)  
     51   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surfr        !: closed sea target rnf surface areas  
     52   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surfe        !: closed sea target empmr surface areas  
    5053 
    5154   !! * Substitutions 
     
    5861CONTAINS 
    5962 
    60    SUBROUTINE dom_clo( cd_cfg, kcfg ) 
     63   SUBROUTINE dom_clo( k_bot ) 
    6164      !!--------------------------------------------------------------------- 
    6265      !!                  ***  ROUTINE dom_clo  *** 
     
    6770      !!                just the thermodynamic processes are applied. 
    6871      !! 
    69       !! ** Action  :   ncsi1(), ncsj1() : south-west closed sea limits (i,j) 
    70       !!                ncsi2(), ncsj2() : north-east Closed sea limits (i,j) 
    71       !!                ncsir(), ncsjr() : Location of runoff 
    72       !!                ncsnr            : number of point where run-off pours 
    73       !!                ncstt            : Type of closed sea 
    74       !!                                   =0 spread over the world ocean 
    75       !!                                   =2 put at location runoff 
    76       !!---------------------------------------------------------------------- 
    77       CHARACTER(len=*), INTENT(in   ) ::   cd_cfg   ! configuration name 
    78       INTEGER         , INTENT(in   ) ::   kcfg     ! configuration identifier  
    79       ! 
    80       INTEGER ::   jc      ! dummy loop indices 
    81       INTEGER ::   isrow   ! local index 
     72      !! ** Action  :   Read closea mask fields from closea_masks.nc and infer 
     73      !!                number of closed seas from closea_mask field. 
     74      !!                closea_mask       : integer values defining closed seas (or groups of closed seas) 
     75      !!                closea_mask_rnf   : integer values defining mappings from closed seas or groups of 
     76      !!                                    closed seas to a runoff area for downwards flux only. 
     77      !!                closea_mask_empmr : integer values defining mappings from closed seas or groups of 
     78      !!                                    closed seas to a runoff area for net fluxes. 
     79      !! 
     80      !!                Python code to generate the closea_masks.nc file from the old-style indices 
     81      !!                definitions is available at TOOLS/DOMAINcfg/make_closea_masks.py 
     82      !!---------------------------------------------------------------------- 
     83      INTEGER, DIMENSION(:,:), INTENT(in)  ::   k_bot   ! ocean last level index (for masking input fields)  
     84      INTEGER ::   inum    ! input file identifier 
     85      INTEGER ::   ierr    ! error code 
     86      INTEGER ::   id      ! netcdf variable ID 
     87      REAL(wp), POINTER, DIMENSION(:,:) :: zdata_in ! temporary real array for input 
    8288      !!---------------------------------------------------------------------- 
    8389      ! 
     
    8692      IF(lwp) WRITE(numout,*)'~~~~~~~' 
    8793      ! 
    88       ! initial values 
    89       ncsnr(:) = 1  ;  ncsi1(:) = 1  ;  ncsi2(:) = 1  ;  ncsir(:,:) = 1 
    90       ncstt(:) = 0  ;  ncsj1(:) = 1  ;  ncsj2(:) = 1  ;  ncsjr(:,:) = 1 
    91       ! 
    92       ! set the closed seas (in data domain indices) 
    93       ! ------------------- 
    94       ! 
    95       IF( cd_cfg == "orca" ) THEN      !==  ORCA configuration  ==! 
    96          ! 
    97          SELECT CASE ( kcfg ) 
    98          !                                           ! ======================= 
    99          CASE ( 1 )                                  !  ORCA_R1 configuration 
    100             !                                        ! ======================= 
    101             IF(lwp) WRITE(numout,*)'   ORCA_R1 closed seas :  only the Caspian Sea' 
    102             ! This dirty section will be suppressed by simplification process: 
    103             ! all this will come back in input files 
    104             ! Currently these hard-wired indices relate to configuration with 
    105             ! extend grid (jpjglo=332) 
    106             isrow = 332 - jpjglo 
    107             ! 
    108             ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea  (spread over the globe) 
    109             ncsi1(1)   = 332  ; ncsj1(1)   = 243 - isrow 
    110             ncsi2(1)   = 344  ; ncsj2(1)   = 275 - isrow 
    111             ncsir(1,1) = 1    ; ncsjr(1,1) = 1 
    112             !                                         
    113             !                                        ! ======================= 
    114          CASE ( 2 )                                  !  ORCA_R2 configuration 
    115             !                                        ! ======================= 
    116             IF(lwp) WRITE(numout,*)'   ORCA_R2 closed seas and lakes : ' 
    117             !                                            ! Caspian Sea 
    118             IF(lwp) WRITE(numout,*)'      Caspian Sea  ' 
    119             ncsnr(1)   =   1  ;  ncstt(1)   =   0           ! spread over the globe 
    120             ncsi1(1)   =  11  ;  ncsj1(1)   = 103 
    121             ncsi2(1)   =  17  ;  ncsj2(1)   = 112 
    122             ncsir(1,1) =   1  ;  ncsjr(1,1) =   1  
    123             !                                            ! Great North American Lakes 
    124             IF(lwp) WRITE(numout,*)'      Great North American Lakes  ' 
    125             ncsnr(2)   =   1  ;  ncstt(2)   =   2           ! put at St Laurent mouth 
    126             ncsi1(2)   =  97  ;  ncsj1(2)   = 107 
    127             ncsi2(2)   = 103  ;  ncsj2(2)   = 111 
    128             ncsir(2,1) = 110  ;  ncsjr(2,1) = 111            
    129             !                                            ! Black Sea (crossed by the cyclic boundary condition) 
    130             IF(lwp) WRITE(numout,*)'      Black Sea  ' 
    131             ncsnr(3:4) =   4  ;  ncstt(3:4) =   2           ! put in Med Sea (north of Aegean Sea) 
    132             ncsir(3:4,1) = 171;  ncsjr(3:4,1) = 106         ! 
    133             ncsir(3:4,2) = 170;  ncsjr(3:4,2) = 106  
    134             ncsir(3:4,3) = 171;  ncsjr(3:4,3) = 105  
    135             ncsir(3:4,4) = 170;  ncsjr(3:4,4) = 105  
    136             ncsi1(3)   = 174  ;  ncsj1(3)   = 107           ! 1 : west part of the Black Sea       
    137             ncsi2(3)   = 181  ;  ncsj2(3)   = 112           !            (ie west of the cyclic b.c.) 
    138             ncsi1(4)   =   2  ;  ncsj1(4)   = 107           ! 2 : east part of the Black Sea  
    139             ncsi2(4)   =   6  ;  ncsj2(4)   = 112           !           (ie east of the cyclic b.c.) 
    140             ! 
    141             !                                        ! ========================= 
    142          CASE ( 025 )                                !  ORCA_R025 configuration 
    143             !                                        ! ========================= 
    144             IF(lwp) WRITE(numout,*)'   ORCA_R025 closed seas : ' 
    145             !                                            ! Caspian Sea 
    146             IF(lwp) WRITE(numout,*)'      Caspian Sea  ' 
    147             ncsnr(1)   = 1    ; ncstt(1)   = 0               ! Caspian + Aral sea 
    148             ncsi1(1)   = 1330 ; ncsj1(1)   = 645 
    149             ncsi2(1)   = 1400 ; ncsj2(1)   = 795 
    150             ncsir(1,1) = 1    ; ncsjr(1,1) = 1 
    151             !                                         
    152             IF(lwp) WRITE(numout,*)'      Azov Sea  ' 
    153             ncsnr(2)   = 1    ; ncstt(2)   = 0               ! Azov Sea  
    154             ncsi1(2)   = 1284 ; ncsj1(2)   = 722 
    155             ncsi2(2)   = 1304 ; ncsj2(2)   = 747 
    156             ncsir(2,1) = 1    ; ncsjr(2,1) = 1 
    157             ! 
    158          END SELECT 
    159          ! 
    160       ELSE                             !==  No closed sea in the configuration  ==! 
    161          ! 
    162          IF(lwp) WRITE(numout,*)'   No closed seas or lakes in the configuration ' 
    163          ! 
     94      CALL wrk_alloc( jpi,jpj, zdata_in) 
     95      !  
     96      ! read the closed seas masks 
     97      ! -------------------------- 
     98      ! 
     99      ALLOCATE( closea_mask(jpi,jpj) , STAT=ierr ) 
     100      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask array') 
     101      ALLOCATE( closea_mask_rnf(jpi,jpj) , STAT=ierr ) 
     102      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_rnf array') 
     103      ALLOCATE( closea_mask_empmr(jpi,jpj) , STAT=ierr ) 
     104      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_empmr array') 
     105 
     106      closea_mask(:,:) = 0 
     107      closea_mask_rnf(:,:) = 0 
     108      closea_mask_empmr(:,:) = 0 
     109 
     110      CALL iom_open( 'closea_masks.nc', inum ) 
     111 
     112      zdata_in(:,:) = 0.0 
     113      CALL iom_get ( inum, jpdom_data, 'closea_mask', zdata_in ) 
     114      closea_mask(:,:) = NINT(zdata_in(:,:)) 
     115      ! necessary because fill values in input fields can confuse things 
     116      ! we can't multiply by tmask here because it isn't defined yet 
     117      WHERE( k_bot(:,:) == 0 ) closea_mask(:,:) = 0 
     118 
     119      id = iom_varid(inum, 'closea_mask_rnf', ldstop = .false.) 
     120      IF(lwp) WRITE(numout,*) 'closea_mask_rnf, id : ',id 
     121      IF( id > 0 ) THEN  
     122         CALL iom_get ( inum, jpdom_data, 'closea_mask_rnf', zdata_in ) 
     123         closea_mask_rnf(:,:) = NINT(zdata_in(:,:)) 
     124         WHERE( k_bot(:,:) == 0 ) closea_mask_rnf(:,:) = 0 
    164125      ENDIF 
    165  
    166       ! convert the position in local domain indices 
    167       ! -------------------------------------------- 
    168       DO jc = 1, jpncs 
    169          ncsi1(jc)   = mi0( ncsi1(jc) ) 
    170          ncsj1(jc)   = mj0( ncsj1(jc) ) 
    171          ! 
    172          ncsi2(jc)   = mi1( ncsi2(jc) )    
    173          ncsj2(jc)   = mj1( ncsj2(jc) )   
    174       END DO 
     126  
     127      id = iom_varid(inum, 'closea_mask_empmr', ldstop = .false.) 
     128      IF( id > 0 ) THEN  
     129         CALL iom_get ( inum, jpdom_data, 'closea_mask_empmr', zdata_in ) 
     130         closea_mask_empmr(:,:) = NINT(zdata_in(:,:)) 
     131         WHERE( k_bot(:,:) == 0 ) closea_mask_empmr(:,:) = 0 
     132      ENDIF 
     133 
     134      CALL iom_close( inum ) 
     135 
     136      ! number of closed seas = global maximum value in closea_mask field 
     137      jncs = maxval(closea_mask(:,:)) 
     138      IF( lk_mpp ) CALL mpp_max(jncs) 
     139      IF( lwp ) WRITE(numout,*) 'Number of closed seas : ',jncs 
     140      ! 
     141      ! number of closed seas rnf mappings = global maximum in closea_mask_rnf field 
     142      jncsr = maxval(closea_mask_rnf(:,:)) 
     143      IF( lk_mpp ) CALL mpp_max(jncsr) 
     144      IF( lwp ) WRITE(numout,*) 'Number of closed seas rnf mappings : ',jncsr 
     145      ! 
     146      ! number of closed seas empmr mappings = global maximum value in closea_mask_empmr field 
     147      jncse = maxval(closea_mask_empmr(:,:)) 
     148      IF( lk_mpp ) CALL mpp_max(jncse) 
     149      IF( lwp ) WRITE(numout,*) 'Number of closed seas empmr mappings : ',jncse 
     150      ! 
     151      CALL wrk_dealloc(jpi,jpj, zdata_in) 
    175152      ! 
    176153   END SUBROUTINE dom_clo 
    177154 
    178155 
    179    SUBROUTINE sbc_clo( kt, cd_cfg, kcfg ) 
     156   SUBROUTINE sbc_clo( kt ) 
    180157      !!--------------------------------------------------------------------- 
    181158      !!                  ***  ROUTINE sbc_clo  *** 
     
    190167      !!---------------------------------------------------------------------- 
    191168      INTEGER         , INTENT(in   ) ::   kt       ! ocean model time step 
    192       CHARACTER(len=*), INTENT(in   ) ::   cd_cfg   ! configuration name 
    193       INTEGER         , INTENT(in   ) ::   kcfg     ! configuration identifier  
    194       ! 
    195       INTEGER             ::   ji, jj, jc, jn   ! dummy loop indices 
     169      ! 
     170      INTEGER             ::   ierr 
     171      INTEGER             ::   jc, jcr, jce   ! dummy loop indices 
    196172      REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon 
    197       REAL(wp)            ::   zze2, ztmp, zcorr     !  
    198       REAL(wp)            ::   zcoef, zcoef1         !  
    199       COMPLEX(wp)         ::   ctmp  
    200       REAL(wp), DIMENSION(jpncs) ::   zfwf   ! 1D workspace 
     173      REAL(wp)            ::   zfwf_total, zcoef, zcoef1         !  
     174      REAL(wp), POINTER, DIMENSION(:)   ::   zfwf, zfwfr, zfwfe     ! freshwater fluxes over closed seas 
     175      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmp2d   ! 2D workspace 
    201176      !!---------------------------------------------------------------------- 
    202177      ! 
    203178      IF( nn_timing == 1 )  CALL timing_start('sbc_clo') 
    204       !                                                   !------------------! 
     179      ! 
     180      CALL wrk_alloc( jncs, zfwf ) 
     181      ! jncsr and jncse can be zero so add 1 to avoid allocating zero-length arrays 
     182      CALL wrk_alloc( jncsr+1, zfwfr ) 
     183      CALL wrk_alloc( jncse+1, zfwfe ) 
     184      CALL wrk_alloc( jpi,jpj, ztmp2d ) 
     185      !                                                   !------------------!  
    205186      IF( kt == nit000 ) THEN                             !  Initialisation  ! 
    206187         !                                                !------------------! 
     
    209190         IF(lwp) WRITE(numout,*)'~~~~~~~' 
    210191 
    211          surf(:) = 0._wp 
    212          ! 
    213          surf(jpncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean 
    214          ! 
    215          !                                        ! surface of closed seas  
    216          DO jc = 1, jpncs 
    217             ctmp = CMPLX( 0.e0, 0.e0, wp ) 
    218             DO jj = ncsj1(jc), ncsj2(jc) 
    219                DO ji = ncsi1(jc), ncsi2(jc) 
    220                   ztmp = e1e2t(ji,jj) * tmask_i(ji,jj) 
    221                   CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
    222                END DO  
    223             END DO  
    224             IF( lk_mpp )   CALL mpp_sum( ctmp ) 
    225             surf(jc) = REAL(ctmp,wp) 
     192         ALLOCATE( surf(jncs+1) , STAT=ierr ) 
     193         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surf array') 
     194         surf(:) = 0.e0_wp 
     195         ! 
     196         ! jncsr can be zero so add 1 to avoid allocating zero-length array 
     197         ALLOCATE( surfr(jncsr+1) , STAT=ierr ) 
     198         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surfr array') 
     199         surfr(:) = 0.e0_wp 
     200         ! 
     201         ! jncse can be zero so add 1 to avoid allocating zero-length array 
     202         ALLOCATE( surfe(jncse+1) , STAT=ierr ) 
     203         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surfe array') 
     204         surfe(:) = 0.e0_wp 
     205         ! 
     206         surf(jncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean 
     207         ! 
     208         !                                        ! surface areas of closed seas  
     209         DO jc = 1, jncs 
     210            ztmp2d(:,:) = 0.e0_wp 
     211            WHERE( closea_mask(:,:) == jc ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:) 
     212            surf(jc) = glob_sum( ztmp2d(:,:) ) 
    226213         END DO 
    227  
    228          IF(lwp) WRITE(numout,*)'     Closed sea surfaces' 
    229          DO jc = 1, jpncs 
    230             IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc) 
     214         ! 
     215         ! jncs+1 : surface area of global ocean, closed seas excluded 
     216         surf(jncs+1) = surf(jncs+1) - SUM(surf(1:jncs)) 
     217         ! 
     218         !                                        ! surface areas of rnf target areas 
     219         DO jcr = 1, jncsr 
     220            ztmp2d(:,:) = 0.e0_wp 
     221            WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:) 
     222            surfr(jcr) = glob_sum( ztmp2d(:,:) ) 
    231223         END DO 
    232  
    233          ! jpncs+1 : surface of sea, closed seas excluded 
    234          DO jc = 1, jpncs 
    235             surf(jpncs+1) = surf(jpncs+1) - surf(jc) 
    236          END DO            
    237          ! 
     224         ! 
     225         !                                        ! surface areas of empmr target areas 
     226         DO jce = 1, jncse 
     227            ztmp2d(:,:) = 0.e0_wp 
     228            WHERE( closea_mask_empmr(:,:) == jcr .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:) 
     229            surfe(jcr) = glob_sum( ztmp2d(:,:) ) 
     230         END DO 
     231         ! 
     232         IF(lwp) WRITE(numout,*)'     Closed sea surface areas (km2)' 
     233         DO jc = 1, jncs 
     234            IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jc, surf(jc) * 1.0e-6 
     235         END DO 
     236         IF(lwp) WRITE(numout,FMT='(A,ES12.2)') 'Global surface area excluding closed seas (km2): ', surf(jncs+1) * 1.0e-6 
     237         ! 
     238         IF(jncsr > 0) THEN 
     239            IF(lwp) WRITE(numout,*)'     Closed sea target rnf surface areas (km2)' 
     240            DO jcr = 1, jncsr 
     241               IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jcr, surfr(jcr) * 1.0e-6 
     242            END DO 
     243         ENDIF 
     244         ! 
     245         IF(jncse > 0) THEN 
     246            IF(lwp) WRITE(numout,*)'     Closed sea target empmr surface areas (km2)' 
     247            DO jce = 1, jncse 
     248               IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jce, surfe(jce) * 1.0e-6 
     249            END DO 
     250         ENDIF 
    238251      ENDIF 
    239       !                                                   !--------------------! 
    240       !                                                   !  update emp        ! 
    241       zfwf = 0.e0_wp                                      !--------------------! 
    242       DO jc = 1, jpncs 
    243          ctmp = CMPLX( 0.e0, 0.e0, wp ) 
    244          DO jj = ncsj1(jc), ncsj2(jc) 
    245             DO ji = ncsi1(jc), ncsi2(jc) 
    246                ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 
    247                CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
    248             END DO   
    249          END DO  
    250          IF( lk_mpp )   CALL mpp_sum( ctmp ) 
    251          zfwf(jc) = REAL(ctmp,wp) 
     252      ! 
     253      !                                                      !--------------------! 
     254      !                                                      !  update emp        ! 
     255      !                                                      !--------------------! 
     256 
     257      zfwf_total = 0._wp 
     258 
     259      ! 
     260      ! 1. Work out total freshwater fluxes over closed seas from EMP - RNF. 
     261      ! 
     262      zfwf(:) = 0.e0_wp            
     263      DO jc = 1, jncs 
     264         ztmp2d(:,:) = 0.e0_wp 
     265         WHERE( closea_mask(:,:) == jc ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) 
     266         zfwf(jc) = glob_sum( ztmp2d(:,:) ) 
    252267      END DO 
    253  
    254       IF( cd_cfg == "orca" .AND. kcfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration 
    255          zze2    = ( zfwf(3) + zfwf(4) ) * 0.5_wp 
    256          zfwf(3) = zze2 
    257          zfwf(4) = zze2 
     268      zfwf_total = SUM(zfwf) 
     269 
     270      ! 
     271      ! 2. Work out total FW fluxes over rnf source areas and add to rnf target areas. If jncsr is zero does not loop. 
     272      !    Where zfwf is negative add flux at specified runoff points and subtract from fluxes for global redistribution. 
     273      !    Where positive leave in global redistribution total. 
     274      ! 
     275      zfwfr(:) = 0.e0_wp            
     276      DO jcr = 1, jncsr 
     277         ! 
     278         ztmp2d(:,:) = 0.e0_wp 
     279         WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) > 0 ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) 
     280         zfwfr(jcr) = glob_sum( ztmp2d(:,:) ) 
     281         IF(lwp) WRITE(numout,*) 'closea runoff output: jcr, zfwfr(jcr), ABS(zfwfr(jcr)/surf(jncs+1) : ',jcr,zfwfr(jcr),ABS(zfwfr(jcr)/surf(jncs+1)) 
     282         ! 
     283         ! The following if avoids the redistribution of the round off 
     284         IF ( ABS(zfwfr(jcr) / surf(jncs+1) ) > rsmall) THEN 
     285            ! 
     286            ! Add residuals to target runoff points if negative and subtract from total to be added globally 
     287            IF(lwp) WRITE(numout,*) 'closea runoff output: passed roundoff test!' 
     288            IF( zfwfr(jcr) < 0.0 ) THEN  
     289               zfwf_total = zfwf_total - zfwfr(jcr) 
     290               zcoef    = zfwfr(jcr) / surfr(jcr) 
     291               zcoef1   = rcp * zcoef 
     292               IF(lwp) WRITE(numout,*) 'closea runoff output: jcr, zcoef: ',jcr, zcoef 
     293               WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0.0) 
     294                  emp(:,:) = emp(:,:) - zcoef 
     295                  qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:) 
     296               ENDWHERE 
     297            ENDIF 
     298            ! 
     299         ENDIF 
     300      END DO 
     301     
     302      ! 
     303      ! 3. Work out total fluxes over empmr source areas and add to empmr target areas. If jncse is zero does not loop.  
     304      ! 
     305      zfwfe(:) = 0.e0_wp            
     306      DO jce = 1, jncse 
     307         ! 
     308         ztmp2d(:,:) = 0.e0_wp 
     309         WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) > 0 ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) 
     310         zfwfe(jce) = glob_sum( ztmp2d(:,:) ) 
     311         ! 
     312         ! The following if avoids the redistribution of the round off 
     313         IF ( ABS(zfwfe(jce) / surf(jncs+1) ) > rsmall) THEN 
     314            ! 
     315            ! Add residuals to runoff points and subtract from total to be added globally 
     316            zfwf_total = zfwf_total - zfwfe(jce) 
     317            zcoef    = zfwfe(jce) / surfe(jce) 
     318            zcoef1   = rcp * zcoef 
     319            WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) == 0.0) 
     320               emp(:,:) = emp(:,:) - zcoef 
     321               qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:) 
     322            ENDWHERE 
     323            ! 
     324         ENDIF 
     325      END DO 
     326     
     327      ! 
     328      ! 4. Spread residual flux over global ocean.  
     329      ! 
     330      ! The following if avoids the redistribution of the round off 
     331      IF ( ABS(zfwf_total / surf(jncs+1) ) > rsmall) THEN 
     332         zcoef    = zfwf_total / surf(jncs+1) 
     333         zcoef1   = rcp * zcoef 
     334         IF(lwp) WRITE(numout,*) 'closea global addition: zfwf_total, zcoef: ', zfwf_total, zcoef 
     335         WHERE( closea_mask(:,:) == 0 ) 
     336            emp(:,:) = emp(:,:) + zcoef 
     337            qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 
     338         ENDWHERE 
    258339      ENDIF 
    259340 
    260       zcorr = 0._wp 
    261  
    262       DO jc = 1, jpncs 
    263          ! 
     341      ! 
     342      ! 5. Subtract area means from emp (and qns) over closed seas to give zero mean FW flux over each sea. 
     343      ! 
     344      DO jc = 1, jncs 
    264345         ! The following if avoids the redistribution of the round off 
    265          IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN 
    266             ! 
    267             IF( ncstt(jc) == 0 ) THEN           ! water/evap excess is shared by all open ocean 
    268                zcoef    = zfwf(jc) / surf(jpncs+1) 
    269                zcoef1   = rcp * zcoef 
    270                emp(:,:) = emp(:,:) + zcoef 
    271                qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 
    272                ! accumulate closed seas correction 
    273                zcorr    = zcorr    + zcoef 
    274                ! 
    275             ELSEIF( ncstt(jc) == 1 ) THEN       ! Excess water in open sea, at outflow location, excess evap shared 
    276                IF ( zfwf(jc) <= 0.e0_wp ) THEN  
    277                    DO jn = 1, ncsnr(jc) 
    278                      ji = mi0(ncsir(jc,jn)) 
    279                      jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 
    280                      IF (      ji > 1 .AND. ji < jpi   & 
    281                          .AND. jj > 1 .AND. jj < jpj ) THEN  
    282                          zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 
    283                          zcoef1     = rcp * zcoef 
    284                          emp(ji,jj) = emp(ji,jj) + zcoef 
    285                          qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj) 
    286                      ENDIF  
    287                    END DO  
    288                ELSE  
    289                    zcoef    = zfwf(jc) / surf(jpncs+1) 
    290                    zcoef1   = rcp * zcoef 
    291                    emp(:,:) = emp(:,:) + zcoef 
    292                    qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 
    293                    ! accumulate closed seas correction 
    294                    zcorr    = zcorr    + zcoef 
    295                ENDIF 
    296             ELSEIF( ncstt(jc) == 2 ) THEN       ! Excess e-p-r (either sign) goes to open ocean, at outflow location 
    297                DO jn = 1, ncsnr(jc) 
    298                   ji = mi0(ncsir(jc,jn)) 
    299                   jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 
    300                   IF(      ji > 1 .AND. ji < jpi    & 
    301                      .AND. jj > 1 .AND. jj < jpj ) THEN  
    302                      zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) ) 
    303                      zcoef1     = rcp * zcoef 
    304                      emp(ji,jj) = emp(ji,jj) + zcoef 
    305                      qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj) 
    306                   ENDIF  
    307                END DO  
    308             ENDIF  
    309             ! 
    310             DO jj = ncsj1(jc), ncsj2(jc) 
    311                DO ji = ncsi1(jc), ncsi2(jc) 
    312                   zcoef      = zfwf(jc) / surf(jc) 
    313                   zcoef1     = rcp * zcoef 
    314                   emp(ji,jj) = emp(ji,jj) - zcoef 
    315                   qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj) 
    316                END DO   
    317             END DO  
    318             ! 
    319          END IF 
    320       END DO  
    321  
    322       IF ( ABS(zcorr) > rsmall ) THEN      ! remove the global correction from the closed seas 
    323          DO jc = 1, jpncs                  ! only if it is large enough 
    324             DO jj = ncsj1(jc), ncsj2(jc) 
    325                DO ji = ncsi1(jc), ncsi2(jc) 
    326                   emp(ji,jj) = emp(ji,jj) - zcorr 
    327                   qns(ji,jj) = qns(ji,jj) + rcp * zcorr * sst_m(ji,jj) 
    328                END DO   
    329              END DO  
    330           END DO 
    331       ENDIF 
     346         IF ( ABS(zfwf(jc) / surf(jncs+1) ) > rsmall) THEN 
     347            ! 
     348            ! Subtract residuals from fluxes over closed sea 
     349            zcoef    = zfwf(jc) / surf(jc) 
     350            zcoef1   = rcp * zcoef 
     351            IF(lwp) WRITE(numout,*) 'closea subtract mean: jc, zfwf(jc), zcoef: ',jc, zfwf(jc), zcoef 
     352            WHERE( closea_mask(:,:) == jc ) 
     353               emp(:,:) = emp(:,:) - zcoef 
     354               qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:) 
     355            ENDWHERE 
     356            ! 
     357         ENDIF 
     358      END DO 
    332359      ! 
    333360      emp (:,:) = emp (:,:) * tmask(:,:,1) 
     
    335362      CALL lbc_lnk( emp , 'T', 1._wp ) 
    336363      ! 
     364      CALL wrk_dealloc( jncs   , zfwf ) 
     365      CALL wrk_dealloc( jncsr+1  , zfwfr ) 
     366      CALL wrk_dealloc( jncse+1  , zfwfe ) 
     367      CALL wrk_dealloc( jpi,jpj, ztmp2d ) 
     368      ! 
    337369      IF( nn_timing == 1 )  CALL timing_stop('sbc_clo') 
    338370      ! 
    339371   END SUBROUTINE sbc_clo 
    340  
    341372 
    342373   SUBROUTINE clo_rnf( p_rnfmsk ) 
     
    353384      !!---------------------------------------------------------------------- 
    354385      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array) 
    355       ! 
    356       INTEGER  ::   jc, jn, ji, jj      ! dummy loop indices 
    357       !!---------------------------------------------------------------------- 
    358       ! 
    359       DO jc = 1, jpncs 
    360          IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows 
    361              DO jn = 1, 4 
    362                 DO jj =    mj0( ncsjr(jc,jn) ), mj1( ncsjr(jc,jn) ) 
    363                    DO ji = mi0( ncsir(jc,jn) ), mi1( ncsir(jc,jn) ) 
    364                       p_rnfmsk(ji,jj) = MAX( p_rnfmsk(ji,jj), 1.0_wp ) 
    365                    END DO 
    366                 END DO 
    367             END DO  
    368          ENDIF  
    369       END DO  
     386      !!---------------------------------------------------------------------- 
     387      ! 
     388      WHERE( ( closea_mask_rnf(:,:) > 0 .or. closea_mask_empmr(:,:) > 0 ) .and. closea_mask(:,:) == 0 ) 
     389         p_rnfmsk(:,:) = MAX( p_rnfmsk(:,:), 1.0_wp ) 
     390      ENDWHERE 
    370391      ! 
    371392   END SUBROUTINE clo_rnf 
     
    383404      !!---------------------------------------------------------------------- 
    384405      INTEGER, DIMENSION(:,:), INTENT(inout) ::   k_top, k_bot   ! ocean first and last level indices 
    385       ! 
    386       INTEGER  ::   jc, ji, jj      ! dummy loop indices 
    387       !!---------------------------------------------------------------------- 
    388       ! 
    389       DO jc = 1, jpncs 
    390          DO jj = ncsj1(jc), ncsj2(jc) 
    391             DO ji = ncsi1(jc), ncsi2(jc) 
    392                k_top(ji,jj) = 0    
    393                k_bot(ji,jj) = 0    
    394             END DO  
    395          END DO  
    396        END DO  
    397        ! 
     406      !!---------------------------------------------------------------------- 
     407      ! 
     408      WHERE( closea_mask(:,:) > 0 ) 
     409         k_top(:,:) = 0    
     410         k_bot(:,:) = 0    
     411      ENDWHERE 
     412      ! 
    398413   END SUBROUTINE clo_bat 
    399414 
    400415   !!====================================================================== 
    401 END MODULE usrdef_closea 
    402  
     416END MODULE closea 
     417 
  • branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r8995 r9017  
    3131   !                                   !!* Namelist namdom : time & space domain * 
    3232   LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time 
     33   LOGICAL , PUBLIC ::   ln_closea      !: =T apply special treatment to closed seas depending on value of nn_closea. 
    3334   INTEGER , PUBLIC ::   nn_closea      !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 
    3435   INTEGER , PUBLIC ::   nn_msh         !: >0  create a mesh-mask file (mesh_mask.nc) 
  • branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r8995 r9017  
    3030   USE trc_oce        ! shared ocean & passive tracers variab 
    3131   USE phycst         ! physical constants 
    32    USE usrdef_closea  ! closed seas 
     32   USE closea         ! closed seas 
    3333   USE domhgr         ! domain: set the horizontal mesh 
    3434   USE domzgr         ! domain: set the vertical mesh 
     
    128128      CALL dom_glo                     ! global domain versus local domain 
    129129      CALL dom_nam                     ! read namelist ( namrun, namdom ) 
    130       CALL dom_clo( cn_cfg, nn_cfg )   ! Closed seas and lake 
    131130      CALL dom_hgr                     ! Horizontal mesh 
    132131      CALL dom_zgr( ik_top, ik_bot )   ! Vertical mesh and bathymetry 
    133       IF( nn_closea == 0 )   CALL clo_bat( ik_top, ik_bot )    !==  remove closed seas or lakes  ==! 
     132      IF( ln_closea ) THEN     
     133          CALL dom_clo( ik_bot )  ! Closed seas and lake  
     134          IF( nn_closea == 0 )   CALL clo_bat( ik_top, ik_bot )    !==  remove closed seas or lakes  ==! 
     135      ENDIF 
    134136      CALL dom_msk( ik_top, ik_bot )   ! Masks 
    135137      ! 
     
    286288         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     & 
    287289         &             ln_cfmeta, ln_iscpl 
    288       NAMELIST/namdom/ ln_linssh, nn_closea, nn_msh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs 
     290      NAMELIST/namdom/ ln_linssh, ln_closea, nn_closea, nn_msh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs 
    289291#if defined key_netcdf4 
    290292      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip 
     
    395397         WRITE(numout,*) '   Namelist namdom : space & time domain' 
    396398         WRITE(numout,*) '      linear free surface (=T)              ln_linssh  = ', ln_linssh 
    397          WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea  = ', nn_closea 
     399         WRITE(numout,*) '      special treatment of closed seas      ln_closea  = ', ln_closea 
     400         IF(ln_closea) WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea  = ', nn_closea 
    398401         WRITE(numout,*) '      create mesh/mask file(s)              nn_msh     = ', nn_msh 
    399402         WRITE(numout,*) '           = 0   no file created           ' 
  • branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r8995 r9017  
    4949   USE bdy_oce   , ONLY: ln_bdy 
    5050   USE usrdef_sbc     ! user defined: surface boundary condition 
    51    USE usrdef_closea  ! user defined: closed sea 
     51   USE closea         ! closed sea 
    5252   ! 
    5353   USE prtctl         ! Print control                    (prt_ctl routine) 
     
    153153         WRITE(numout,*) '         runoff / runoff mouths                     ln_rnf        = ', ln_rnf 
    154154         WRITE(numout,*) '         iceshelf formulation                       ln_isf        = ', ln_isf 
    155          WRITE(numout,*) '         closed sea (=0/1) (set in namdom)          nn_closea     = ', nn_closea 
     155         WRITE(numout,*) '         special treatment of closed seas (set in namdom)  ln_closea     = ', ln_closea 
     156         IF(ln_closea) WRITE(numout,*) '         closed sea treatment (=0/1) (set in namdom) nn_closea     = ', nn_closea 
    156157         WRITE(numout,*) '         nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm 
    157158         WRITE(numout,*) '         surface wave                               ln_wave       = ', ln_wave 
     
    446447      ! treatment of closed sea in the model domain   (update freshwater fluxes) 
    447448      ! Should not be ran if ln_diurnal_only 
    448       IF( .NOT.ln_diurnal_only .AND. nn_closea == 1 )   CALL sbc_clo( kt, cn_cfg, nn_cfg )    
     449      IF( .NOT.ln_diurnal_only .AND. ln_closea .AND. nn_closea == 1 )   CALL sbc_clo( kt )    
    449450 
    450451!RBbug do not understand why see ticket 667 
  • branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r8995 r9017  
    2121   USE sbcisf         ! PM we could remove it I think 
    2222   USE eosbn2         ! Equation Of State 
    23    USE usrdef_closea  ! closed seas 
     23   USE closea         ! closed seas 
    2424   ! 
    2525   USE in_out_manager ! I/O manager 
     
    513513      CALL iom_close( inum )                                      ! close file 
    514514      ! 
    515       IF( nn_closea == 1 )   CALL clo_rnf( rnfmsk )               ! closed sea inflow set as ruver mouth 
     515      IF( ln_closea .and. nn_closea == 1 )   CALL clo_rnf( rnfmsk )   ! closed sea inflow set as river mouth 
    516516      ! 
    517517      rnfmsk_z(:)   = 0._wp                                       ! vertical structure 
Note: See TracChangeset for help on using the changeset viewer.