Changeset 13024
- Timestamp:
- 2020-06-03T16:26:23+02:00 (5 years ago)
- Location:
- utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg
- Files:
-
- 28 added
- 25 edited
Legend:
- Unmodified
- Added
- Removed
-
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/namelist_cfg
r12414 r13024 15 15 &namdom ! space and time domain (bathymetry, mesh, timestep) 16 16 !----------------------------------------------------------------------- 17 nn_bathy = 1 ! compute analyticaly (=0) or read (=1) the bathymetry file 18 ! or compute (2) from external bathymetry 19 nn_interp = 1 ! type of interpolation (nn_bathy =2) 20 cn_topo = 'bathymetry_ORCA12_V3.3.nc' ! external topo file (nn_bathy =2) 21 cn_bath = 'Bathymetry' ! topo name in file (nn_bathy =2) 22 cn_lon = 'nav_lon' ! lon name in file (nn_bathy =2) 23 cn_lat = 'nav_lat' ! lat name in file (nn_bathy =2) 24 rn_scale = 1 25 rn_bathy = 0. ! value of the bathymetry. if (=0) bottom flat at jpkm1 26 jphgr_msh = 0 ! type of horizontal mesh 27 ppglam0 = 999999.0 ! longitude of first raw and column T-point (jphgr_msh = 1) 28 ppgphi0 = 999999.0 ! latitude of first raw and column T-point (jphgr_msh = 1) 29 ppe1_deg = 999999.0 ! zonal grid-spacing (degrees) 30 ppe2_deg = 999999.0 ! meridional grid-spacing (degrees) 31 ppe1_m = 999999.0 ! zonal grid-spacing (degrees) 32 ppe2_m = 999999.0 ! meridional grid-spacing (degrees) 33 ppsur = -4762.96143546300 ! ORCA r4, r2 and r05 coefficients 34 ppa0 = 255.58049070440 ! (default coefficients) 35 ppa1 = 245.58132232490 ! 36 ppkth = 21.43336197938 ! 37 ppacr = 3.0 ! 38 ppdzmin = 999999. ! Minimum vertical spacing 39 pphmax = 999999. ! Maximum depth 40 ldbletanh = .FALSE. ! Use/do not use double tanf function for vertical coordinates 41 ppa2 = 999999. ! Double tanh function parameters 42 ppkth2 = 999999. ! 43 ppacr2 = 999999. ! 44 / 45 !----------------------------------------------------------------------- 46 &namcrs ! Grid coarsening for dynamics output and/or 47 ! passive tracer coarsened online simulations 48 !----------------------------------------------------------------------- 49 / 50 !----------------------------------------------------------------------- 51 &namtsd ! data : Temperature & Salinity 52 !----------------------------------------------------------------------- 53 / 54 !----------------------------------------------------------------------- 55 &namsbc ! Surface Boundary Condition (surface module) 56 !----------------------------------------------------------------------- 57 / 58 !----------------------------------------------------------------------- 59 &namsbc_core ! namsbc_core CORE bulk formulae 60 !----------------------------------------------------------------------- 61 / 62 !----------------------------------------------------------------------- 63 &namtra_qsr ! penetrative solar radiation 64 !----------------------------------------------------------------------- 65 / 66 !----------------------------------------------------------------------- 67 &namsbc_rnf ! runoffs namelist surface boundary condition 68 !----------------------------------------------------------------------- 69 / 70 !----------------------------------------------------------------------- 71 &namsbc_ssr ! surface boundary condition : sea surface restoring 72 !----------------------------------------------------------------------- 73 / 74 !----------------------------------------------------------------------- 75 &namsbc_alb ! albedo parameters 76 !----------------------------------------------------------------------- 77 / 78 !----------------------------------------------------------------------- 79 &namberg ! iceberg parameters 80 !----------------------------------------------------------------------- 81 / 82 !----------------------------------------------------------------------- 83 &namlbc ! lateral momentum boundary condition 84 !----------------------------------------------------------------------- 85 / 86 !----------------------------------------------------------------------- 87 &nambfr ! bottom friction 88 !----------------------------------------------------------------------- 89 / 90 !----------------------------------------------------------------------- 91 &nambbc ! bottom temperature boundary condition (default: NO) 92 !----------------------------------------------------------------------- 93 ln_trabbc = .true. ! Apply a geothermal heating at the ocean bottom 94 / 95 !----------------------------------------------------------------------- 96 &nambbl ! bottom boundary layer scheme 97 !----------------------------------------------------------------------- 17 98 / 18 99 !----------------------------------------------------------------------- 19 100 &namcfg ! parameters of the configuration 20 101 !----------------------------------------------------------------------- 102 ! 103 ln_e3_dep = .true. ! =T : e3=dk[depth] in discret sens. 104 ! ! ===>>> will become the only possibility in v4.0 105 ! ! =F : e3 analytical derivative of depth function 106 ! ! only there for backward compatibility test with v3.6 107 ! ! 108 cp_cfg = "orca" ! name of the configuration 109 jp_cfg = 2 ! resolution of the configuration 110 jpidta = 182 ! 1st lateral dimension ( >= jpi ) 111 jpjdta = 149 ! 2nd " " ( >= jpj ) 112 jpkdta = 31 ! number of levels ( >= jpk ) 113 jpiglo = 182 ! 1st dimension of global domain --> i =jpidta 114 jpjglo = 149 ! 2nd - - --> j =jpjdta 115 jperio = 4 ! lateral cond. type (between 0 and 6) 116 ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present 117 ! in netcdf input files, as the start j-row for reading 118 ln_domclo = .false. ! computation of closed sea masks (see namclo) 21 119 / 22 120 !----------------------------------------------------------------------- 23 121 &namzgr ! vertical coordinate (default: NO selection) 24 122 !----------------------------------------------------------------------- 123 !----------------------------------------------------------------------- 124 ln_zco = .false. ! z-coordinate - full steps 125 ln_zps = .true. ! z-coordinate - partial steps 126 ln_sco = .false. ! s- or hybrid z-s-coordinate 127 ln_isfcav = .false. ! ice shelf cavity (T: see namzgr_isf) 128 ln_linssh = .true. ! linear free surface 25 129 / 26 130 !----------------------------------------------------------------------- … … 39 143 &namlbc ! lateral momentum boundary condition (default: NO selection) 40 144 !----------------------------------------------------------------------- 145 jpni = 0 146 jpnj=0 41 147 / 42 148 !----------------------------------------------------------------------- 43 149 &namagrif ! AGRIF zoom ("key_agrif") 150 ln_bry_south = .TRUE. 44 151 !----------------------------------------------------------------------- 45 152 / -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/namelist_ref
r12414 r13024 40 40 &namdom ! space and time domain (bathymetry, mesh, timestep) 41 41 !----------------------------------------------------------------------- 42 nn_bathy = 1 ! compute (=0) or read (=1) the bathymetry file 42 nn_bathy = 1 ! compute analyticaly (=0) or read (=1) the bathymetry file 43 ! or compute (2) from external bathymetry 44 nn_interp = 1 ! type of interpolation (nn_bathy =2) 45 cn_topo = 'bathymetry_ORCA12_V3.3.nc' ! external topo file (nn_bathy =2) 46 cn_bath = 'Bathymetry' ! topo name in file (nn_bathy =2) 47 cn_lon = 'nav_lon' ! lon name in file (nn_bathy =2) 48 cn_lat = 'nav_lat' ! lat name in file (nn_bathy =2) 43 49 rn_bathy = 0. ! value of the bathymetry. if (=0) bottom flat at jpkm1 44 50 nn_msh = 0 ! create (=1) a mesh file or not (=0) … … 73 79 ppkth2 = 48.029893720000 ! 74 80 ppacr2 = 13.000000000000 ! 75 / 81 // 76 82 !----------------------------------------------------------------------- 77 83 &namcfg ! parameters of the configuration … … 82 88 ! ! =F : e3 analytical derivative of depth function 83 89 ! ! only there for backward compatibility test with v3.6 84 ! 85 cp_cfg = "default" ! name of the configuration 86 cp_cfz = "no zoom" ! name of the zoom of configuration 87 jp_cfg = 0 ! resolution of the configuration 88 jpkdta = 31 ! number of levels ( >= jpk ) 89 jpiglo = 10 ! 1st dimension of global domain --> i =jpidta 90 jpjglo = 12 ! 2nd - - --> j =jpjdta 91 jperio = 0 ! lateral cond. type (between 0 and 6) 92 ! = 0 closed ; = 1 cyclic East-West 93 ! = 2 equatorial symmetric ; = 3 North fold T-point pivot 94 ! = 4 cyclic East-West AND North fold T-point pivot 95 ! = 5 North fold F-point pivot 96 ! = 6 cyclic East-West AND North fold F-point pivot 90 ! ! 91 cp_cfg = "orca" ! name of the configuration 92 jp_cfg = 2 ! resolution of the configuration 93 jpidta = 182 ! 1st lateral dimension ( >= jpi ) 94 jpjdta = 149 ! 2nd " " ( >= jpj ) 95 jpkdta = 31 ! number of levels ( >= jpk ) 96 jpiglo = 182 ! 1st dimension of global domain --> i =jpidta 97 jpjglo = 149 ! 2nd - - --> j =jpjdta 98 jperio = 4 ! lateral cond. type (between 0 and 6) 97 99 ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present 98 100 ! in netcdf input files, as the start j-row for reading … … 182 184 rn_sponge_tra = 2880. ! coefficient for tracer sponge layer [m2/s] 183 185 rn_sponge_dyn = 2880. ! coefficient for dynamics sponge layer [m2/s] 184 ln_chk_bathy = .false. ! =T check the parent bathymetry 186 ln_chk_bathy = .FALSE. ! 187 npt_connect = 2 188 npt_copy = 2 189 ln_bry_south = .TRUE. 190 / 191 !----------------------------------------------------------------------- 192 &nam_tide ! tide parameters ("key_tide") 193 !----------------------------------------------------------------------- 194 ln_tide_pot = .true. ! use tidal potential forcing 195 ln_tide_ramp= .false. ! 196 rdttideramp = 0. ! 197 clname(1) = 'DUMMY' ! name of constituent - all tidal components must be set in namelist_cfg 185 198 / 186 199 !----------------------------------------------------------------------- … … 210 223 rn_ice_sal = 10. ! si3 only: -- salinity -- 211 224 rn_ice_age = 30. ! si3 only: -- age -- 212 !213 ln_tra_dmp =.false. ! open boudaries conditions for tracers214 ln_dyn3d_dmp =.false. ! open boundary condition for baroclinic velocities215 rn_time_dmp = 1. ! Damping time scale in days216 rn_time_dmp_out = 1. ! Outflow damping time scale217 nn_rimwidth = 10 ! width of the relaxation zone218 ln_vol = .false. ! total volume correction (see nn_volctl parameter)219 nn_volctl = 1 ! = 0, the total water flux across open boundaries is zero220 nb_jpk_bdy = -1 ! number of levels in the bdy data (set < 0 if consistent with planned run)221 /222 !-----------------------------------------------------------------------223 &namnc4 ! netcdf4 chunking and compression settings ("key_netcdf4")224 !-----------------------------------------------------------------------225 nn_nchunks_i = 4 ! number of chunks in i-dimension226 nn_nchunks_j = 4 ! number of chunks in j-dimension227 nn_nchunks_k = 31 ! number of chunks in k-dimension228 ! ! setting nn_nchunks_k = jpk will give a chunk size of 1 in the vertical which229 ! ! is optimal for postprocessing which works exclusively with horizontal slabs230 ln_nc4zip = .true. ! (T) use netcdf4 chunking and compression231 ! ! (F) ignore chunking information and produce netcdf3-compatible files232 225 / 233 226 !----------------------------------------------------------------------- … … 238 231 nn_buffer = 0 ! size in bytes of exported buffer ('B' case), 0 no exportation 239 232 ln_nnogather = .true. ! activate code to avoid mpi_allgather use at the northfold 240 jpni = 0! jpni number of processors following i (set automatically if < 1)241 jpnj = 0! jpnj number of processors following j (set automatically if < 1)233 jpni = 1 ! jpni number of processors following i (set automatically if < 1) 234 jpnj = 1 ! jpnj number of processors following j (set automatically if < 1) 242 235 / 243 236 !----------------------------------------------------------------------- -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_boundary_connections.F90
r12414 r13024 13 13 call Agrif_Bc_variable(e3t_copy_id, procname = connect_e3t_copy) 14 14 15 ! Agrif_UseSpecialValue = .TRUE. 16 ! Agrif_SpecialValue = 0. 17 ! call Agrif_Bc_variable(e3t_connect_id, procname = connect_e3t_connect) 18 ! Agrif_UseSpecialValue = .FALSE. 19 15 Allocate(e3t_interp(jpi,jpj,jpk)) 16 e3t_interp = -10. 17 Agrif_UseSpecialValue = .TRUE. 18 Agrif_SpecialValue = 0. 19 call Agrif_Bc_variable(e3t_connect_id, procname = connect_e3t_connect) 20 Agrif_UseSpecialValue = .FALSE. 21 22 ! Call Agrif_make_connection() 20 23 21 24 Agrif_SpecialValue = 0. 22 25 Agrif_UseSpecialValue = ln_spc_dyn 23 26 ! 24 CALL Agrif_Bc_variable( e3u_id, procname=connect_e3u )25 CALL Agrif_Bc_variable( e3v_id, procname=connect_e3v )27 ! CALL Agrif_Bc_variable( e3u_id, procname=connect_e3u ) 28 ! CALL Agrif_Bc_variable( e3v_id, procname=connect_e3v ) 26 29 ! 27 30 Agrif_UseSpecialValue = .FALSE. … … 83 86 ELSE 84 87 mbkt(i1:i2,j1:j2) = nint(ptab(i1:i2,j1:j2)) 85 86 88 WHERE (mbkt(i1:i2,j1:j2)==0) 87 89 ssmask(i1:i2,j1:j2) = 0. 90 ELSEWHERE 91 ssmask(i1:i2,j1:j2) = 1. 88 92 END WHERE 89 93 … … 136 140 do ji=i1,i2 137 141 bathy_local (ji,jj) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj) 138 print *,'ji = ',ji,jj,bathy_local(ji,jj),ptab(ji,jj,jpk+1) 139 enddo 140 enddo 141 142 ! DO jk=1,jpk 143 ! DO jj=j1,j2 144 ! DO ji=i1,i2 145 ! e3t_0(ji,jj,jk) = MAX(ptab(ji,jj,jk),MIN(e3zps_min, e3t_1d(jk)*e3zps_rat )) 146 ! e3t_0(ji,jj,jk) = MIN(e3t_0(ji,jj,jk),e3t_1d(jk)) 147 ! ENDDO 148 ! ENDDO 149 ! ENDDO 142 enddo 143 enddo 144 145 DO jk=1,jpk 146 DO jj=j1,j2 147 DO ji=i1,i2 148 if (e3t_interp(ji,jj,jk) == -10) then ! the connection has not yet been done 149 e3t_interp(ji,jj,jk) = MAX(ptab(ji,jj,jk),MIN(e3zps_min, e3t_1d(jk)*e3zps_rat )) 150 e3t_interp(ji,jj,jk) = MIN(e3t_interp(ji,jj,jk),e3t_1d(jk)) 151 e3t_0(ji,jj,jk) = ztabramp(ji,jj)*e3t_0(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_interp(ji,jj,jk) 152 endif 153 ENDDO 154 ENDDO 155 ENDDO 150 156 ENDIF 151 157 ! -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_connection.F90
r12414 r13024 7 7 !!---------------------------------------------------------------------- 8 8 INTEGER :: ji, jj, ind1, ind2 9 INTEGER :: ispongearea 9 INTEGER :: ispongearea, istart 10 10 REAL(wp) :: z1_spongearea 11 11 !!---------------------------------------------------------------------- … … 15 15 16 16 ALLOCATE(ztabramp(jpi,jpj)) 17 ispongearea = 1 + npt_connect * Agrif_irhox() 17 ispongearea = 1 + npt_connect * Agrif_irhox() 18 istart = npt_copy * Agrif_irhox() 18 19 z1_spongearea = 1._wp / REAL( ispongearea ) 19 20 … … 21 22 22 23 ! --- West --- ! 23 IF( ( nbondi == -1) .OR. (nbondi == 2)) THEN24 ind1 = 1+nbghostcells 25 ind2 = 1+nbghostcells+ ispongearea24 IF( ((nbondi == -1) .OR. (nbondi == 2) ).AND. .NOT. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6)) THEN 25 ind1 = 1+nbghostcells + istart 26 ind2 = ind1 + ispongearea 26 27 DO jj = 1, jpj 27 28 DO ji = ind1, ind2 … … 32 33 33 34 ! --- East --- ! 34 IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 35 ind1 = nlci - nbghostcells - ispongearea 36 ind2 = nlci - nbghostcells 35 IF( ((nbondi == 1) .OR. (nbondi == 2) ).AND. .NOT. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6)) THEN 36 !ind1 = nlci - nbghostcells - ispongearea 37 ind2 = nlci - nbghostcells - istart 38 ind1 = ind2 -ispongearea 39 40 37 41 DO jj = 1, jpj 38 42 DO ji = ind1, ind2 … … 43 47 44 48 ! --- South --- ! 45 IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 46 ind1 = 1+nbghostcells 47 ind2 = 1+nbghostcells + ispongearea 49 IF(( (nbondj == -1) .OR. (nbondj == 2) ).AND.(ln_bry_south)) THEN 50 ! ind1 = 1+nbghostcells 51 ! ind2 = 1+nbghostcells + ispongearea 52 ind1 = 1+nbghostcells + istart 53 ind2 = ind1 + ispongearea 48 54 DO jj = ind1, ind2 49 55 DO ji = 1, jpi … … 55 61 ! --- North --- ! 56 62 IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 57 ind1 = nlcj - nbghostcells - ispongearea 58 ind2 = nlcj - nbghostcells 63 ! ind1 = nlcj - nbghostcells - ispongearea 64 ! ind2 = nlcj - nbghostcells 65 66 ind2 = nlcj - nbghostcells - istart 67 ind1 = ind2 -ispongearea 68 59 69 DO jj = ind1, ind2 60 70 DO ji = 1, jpi … … 67 77 END SUBROUTINE Agrif_connection 68 78 79 SUBROUTINE Agrif_make_connection 80 use dom_oce 81 use agrif_parameters 82 !!---------------------------------------------------------------------- 83 !! *** ROUTINE Agrif_Sponge *** 84 !!---------------------------------------------------------------------- 85 INTEGER :: ji, jj, ind1, ind2 86 INTEGER :: ispongearea, istart 87 REAL(wp) :: z1_spongearea 88 !!---------------------------------------------------------------------- 89 ! 90 ! Define ramp from boundaries towards domain interior at T-points 91 ! Store it in ztabramp 92 93 ispongearea = 1 + npt_connect * Agrif_irhox() 94 istart = npt_copy * Agrif_irhox() 95 96 ! --- West --- ! 97 IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 98 99 ind1 = 1+nbghostcells + istart 100 ind2 = ind1 + ispongearea 101 DO jk=1,jpk 102 DO jj = 1, jpj 103 DO ji = ind1, ind2 104 105 ! print *,'VAL = ',ztabramp(ji,jj)*e3t_interp(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_0(ji,jj,jk), & 106 ! e3t_0(ji,jj,jk) 107 e3t_0(ji,jj,jk) = ztabramp(ji,jj)*e3t_interp(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_0(ji,jj,jk) 108 ENDDO 109 ENDDO 110 ENDDO 111 ENDIF 112 113 ! --- East --- ! 114 IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 115 ind2 = nlci - nbghostcells - istart 116 ind1 = ind2 -ispongearea 117 DO jk=1,jpk 118 DO jj = 1, jpj 119 DO ji = ind1, ind2 120 e3t_0(ji,jj,jk) = ztabramp(ji,jj)*e3t_interp(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_0(ji,jj,jk) 121 ENDDO 122 ENDDO 123 ENDDO 124 ENDIF 125 126 ! --- South --- ! 127 IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 128 ind1 = 1+nbghostcells + istart 129 ind2 = ind1 + ispongearea 130 DO jk=1,jpk 131 DO jj = ind1, ind2 132 DO ji = 1, jpi 133 e3t_0(ji,jj,jk) = ztabramp(ji,jj)*e3t_interp(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_0(ji,jj,jk) 134 END DO 135 ENDDO 136 ENDDO 137 ENDIF 138 139 ! --- North --- ! 140 IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 141 142 ind2 = nlcj - nbghostcells - istart 143 ind1 = ind2 -ispongearea 144 DO jk=1,jpk 145 DO jj = ind1, ind2 146 DO ji = 1, jpi 147 e3t_0(ji,jj,jk) = ztabramp(ji,jj)*e3t_interp(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_0(ji,jj,jk) 148 END DO 149 ENDDO 150 ENDDO 151 ENDIF 152 ! 153 ! 154 END SUBROUTINE Agrif_make_connection 155 69 156 #else 70 157 subroutine agrif_connection_empty -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_domzgr.F90
r12414 r13024 23 23 ! 24 24 !!---------------------------------------------------------------------- 25 INTEGER :: ji,jj 25 26 ! 26 27 western_side = (nb == 1).AND.(ndir == 1) … … 30 31 IF( before) THEN 31 32 ptab(i1:i2,j1:j2) = bathy(i1:i2,j1:j2) 33 34 do jj=j1,j2 35 do ji=i1,i2 36 ptab(ji,jj) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj) 37 enddo 38 enddo 32 39 ELSE 33 40 bathy(i1:i2,j1:j2)=ptab -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_parameters.F90
r12414 r13024 10 10 INTEGER :: npt_copy 11 11 INTEGER :: npt_connect 12 12 LOGICAL :: ln_bry_south = .TRUE. 13 13 REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: ztabramp 14 REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:,:) :: e3t_interp 14 15 15 16 end module agrif_parameters -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_recompute_scalefactors.f90
r12414 r13024 9 9 ! Scale factors and depth at U-, V-, UW and VW-points 10 10 DO jk = 1, jpk ! initialisation to z-scale factors 11 !e3u_0 (:,:,jk) = e3t_1d(jk)12 !e3v_0 (:,:,jk) = e3t_1d(jk)11 e3u_0 (:,:,jk) = e3t_1d(jk) 12 e3v_0 (:,:,jk) = e3t_1d(jk) 13 13 e3uw_0(:,:,jk) = e3w_1d(jk) 14 14 e3vw_0(:,:,jk) = e3w_1d(jk) … … 18 18 DO jj = 1, jpjm1 19 19 DO ji = 1, jpim1 ! vector opt. 20 !e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )21 !e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )20 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 21 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 22 22 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 23 23 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) … … 27 27 IF ( ln_isfcav ) THEN 28 28 ! (ISF) define e3uw (adapted for 2 cells in the water column) 29 print *,'NOT READY SINCE:' 30 print *,'MBATHY HAS NOT BEEN CORRECTED / UPDATED' 31 print *,'EVEN NOT COMPUTED IN THE CASE ln_read_cfg = .TRUE.' 32 STOP 29 33 DO jj = 2, jpjm1 30 34 DO ji = 2, jpim1 ! vector opt. -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_update.F90
r12414 r13024 6 6 7 7 if (agrif_root()) return 8 8 9 call agrif_update_variable(bottom_level_id,locupdate=(/npt_copy,0/),procname = update_bottom_level) 10 9 11 10 12 Agrif_UseSpecialValueInUpdate = .TRUE. … … 14 16 Agrif_UseSpecialValueInUpdate = .FALSE. 15 17 16 call agrif_update_variable(e3u_id,procname = update_e3u)17 call agrif_update_variable(e3v_id,procname = update_e3v)18 !call agrif_update_variable(e3u_id,procname = update_e3u) 19 !call agrif_update_variable(e3v_id,procname = update_e3v) 18 20 19 21 end subroutine agrif_update_all … … 46 48 WHERE (mbkt(i1:i2,j1:j2)==0) 47 49 ssmask(i1:i2,j1:j2) = 0. 50 ELSEWHERE 51 ssmask(i1:i2,j1:j2) = 1. 48 52 END WHERE 49 53 … … 69 73 DO jj=j1,j2 70 74 DO ji=i1,i2 71 if (mbkt(ji,jj) < jk) then75 if (mbkt(ji,jj) <= jk) then 72 76 tabres(ji,jj,jk) = e3t_0(ji,jj,jk) 73 77 else … … 81 85 DO jj=j1,j2 82 86 DO ji=i1,i2 83 if (mbkt(ji,jj) < jk) then 87 if (mbkt(ji,jj) <= jk) then 88 e3t_0(ji,jj,jk) = MAX(tabres(ji,jj,jk),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 89 else 84 90 e3t_0(ji,jj,jk) = e3t_1d(jk) 85 else86 e3t_0(ji,jj,jk) = MAX(tabres(ji,jj,jk),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat))87 91 endif 88 92 END DO … … 113 117 do jj=j1,j2 114 118 do ji=i1,i2 115 if (min(mbkt(ji,jj),mbkt(ji+1,jj))< jk) then119 if (min(mbkt(ji,jj),mbkt(ji+1,jj))<=jk) then 116 120 tabres(ji,jj,jk) = zrhoy * e2u(ji,jj) * MIN(e3zps_min,e3t_1d(jk)*e3zps_rat) 117 121 else … … 125 129 DO jj=j1,j2 126 130 DO ji=i1,i2 127 if (min(mbkt(ji,jj),mbkt(ji+1,jj))< jk) then128 e3u_0(ji,jj,jk)= e3t_1d(jk)131 if (min(mbkt(ji,jj),mbkt(ji+1,jj))<=jk) then 132 e3u_0(ji,jj,jk)=MAX(tabres(ji,jj,jk) / e2u(ji,jj),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 129 133 else 130 e3u_0(ji,jj,jk) = MAX(tabres(ji,jj,jk) / e2u(ji,jj),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat))134 e3u_0(ji,jj,jk) = e3t_1d(jk) 131 135 endif 132 136 END DO … … 157 161 do jj=j1,j2 158 162 do ji=i1,i2 159 if (min(mbkt(ji,jj),mbkt(ji,jj+1))< jk) then163 if (min(mbkt(ji,jj),mbkt(ji,jj+1))<=jk) then 160 164 tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * MIN(e3zps_min,e3t_1d(jk)*e3zps_rat) 161 165 else … … 169 173 DO jj=j1,j2 170 174 DO ji=i1,i2 171 if (min(mbkt(ji,jj),mbkt(ji,jj+1))< jk) then172 e3v_0(ji,jj,jk)= e3t_1d(jk)175 if (min(mbkt(ji,jj),mbkt(ji,jj+1))<=jk) then 176 e3v_0(ji,jj,jk)=MAX(tabres(ji,jj,jk) / e1v(ji,jj),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 173 177 else 174 e3v_0(ji,jj,jk) = MAX(tabres(ji,jj,jk) / e1v(ji,jj),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat))178 e3v_0(ji,jj,jk) = e3t_1d(jk) 175 179 endif 176 180 END DO -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_user.F90
r12414 r13024 1 1 #if defined key_agrif 2 subroutineagrif_initworkspace()2 SUBROUTINE agrif_initworkspace() 3 3 !!---------------------------------------------------------------------- 4 4 !! *** ROUTINE Agrif_InitWorkspace *** 5 5 !!---------------------------------------------------------------------- 6 end subroutine agrif_initworkspace 7 SUBROUTINE Agrif_InitValues 6 END SUBROUTINE agrif_initworkspace 7 8 SUBROUTINE Agrif_InitValues 8 9 !!---------------------------------------------------------------------- 9 10 !! *** ROUTINE Agrif_InitValues *** … … 11 12 !! ** Purpose :: Declaration of variables to be interpolated 12 13 !!---------------------------------------------------------------------- 13 USE Agrif_Util 14 USE dom_oce 15 USE nemogcm 16 USE domain 17 !! 18 IMPLICIT NONE 19 20 21 CALL nemo_init !* Initializations of each fine grid 22 23 CALL dom_nam 24 CALL cfg_write ! create the configuration file 25 26 END SUBROUTINE Agrif_InitValues 27 28 SUBROUTINE Agrif_InitValues_cont 29 30 use dom_oce 31 integer :: irafx, irafy 32 logical :: ln_perio 33 integer nx,ny 34 35 irafx = agrif_irhox() 36 irafy = agrif_irhoy() 37 38 nx=nlci ; ny=nlcj 14 USE Agrif_Util 15 USE dom_oce 16 USE nemogcm 17 USE domain 18 !! 19 IMPLICIT NONE 20 21 ! No temporal refinement 22 CALL Agrif_Set_coeffreft(1) 23 24 CALL nemo_init !* Initializations of each fine grid 25 26 CALL dom_nam 27 28 END SUBROUTINE Agrif_InitValues 29 30 SUBROUTINE Agrif_InitValues_cont 31 !!---------------------------------------------------------------------- 32 !! *** ROUTINE Agrif_InitValues_cont *** 33 !! 34 !! ** Purpose :: Initialisation of variables to be interpolated 35 !!---------------------------------------------------------------------- 36 USE dom_oce 37 USE lbclnk 38 !! 39 IMPLICIT NONE 40 ! 41 INTEGER :: nx, ny 42 INTEGER :: irafx, irafy 43 LOGICAL :: ln_perio 44 ! 45 irafx = agrif_irhox() 46 irafy = agrif_irhoy() 47 48 nx = nlci ; ny = nlcj 39 49 40 50 ! IF(jperio /=1 .AND. jperio/=4 .AND. jperio/=6 ) THEN … … 44 54 ! nbghostcellsfine_tot_y= nbghostcellsfine_y +1 45 55 ! ELSE 46 ! nx = (nbcellsx)+2*nbghostcellsfine_x 56 ! nx = (nbcellsx)+2*nbghostcellsfine_x 47 57 ! ny = (nbcellsy)+2*nbghostcellsfine+2 48 58 ! nbghostcellsfine_tot_x= 1 … … 53 63 ! nx = nbcellsx+irafx 54 64 ! ny = nbcellsy+irafy 55 56 WRITE(*,*) ' ' 57 WRITE(*,*)'Size of the High resolution grid: ',nx,' x ',ny 58 WRITE(*,*) ' ' 59 60 call agrif_init_lonlat() 61 ln_perio=.FALSE. 62 if( jperio ==1 .OR. jperio==2 .OR. jperio==4) ln_perio=.TRUE. 63 64 where (glamt < -180) glamt = glamt +360. 65 if (ln_perio) then 66 glamt(1,:)=glamt(nx-1,:) 67 glamt(nx,:)=glamt(2,:) 68 endif 69 70 where (glamu < -180) glamu = glamu +360. 71 if (ln_perio) then 72 glamu(1,:)=glamu(nx-1,:) 73 glamu(nx,:)=glamu(2,:) 74 endif 75 76 where (glamv < -180) glamv = glamv +360. 77 if (ln_perio) then 78 glamv(1,:)=glamv(nx-1,:) 79 glamv(nx,:)=glamv(2,:) 80 endif 81 82 where (glamf < -180) glamf = glamf +360. 83 if (ln_perio) then 84 glamf(1,:)=glamf(nx-1,:) 85 glamf(nx,:)=glamf(2,:) 86 endif 87 88 call agrif_init_scales() 89 90 91 END SUBROUTINE Agrif_InitValues_cont 92 93 94 subroutine agrif_declare_var() 95 use par_oce 96 use dom_oce 97 use agrif_profiles 98 use agrif_parameters 99 100 IMPLICIT NONE 101 102 integer :: ind1, ind2, ind3 103 integer nx,ny 104 integer nbghostcellsfine_tot_x, nbghostcellsfine_tot_y 105 INTEGER :: irafx 106 !!---------------------------------------------------------------------- 107 108 ! 1. Declaration of the type of variable which have to be interpolated 109 !--------------------------------------------------------------------- 110 nx=nlci ; ny=nlcj 111 112 !ind2 = nbghostcellsfine_tot_x + 1 113 !ind3 = nbghostcellsfine_tot_y + 1 114 ind2 = 2 + nbghostcells 115 ind3 = ind2 116 nbghostcellsfine_tot_x=nbghostcells+1 117 nbghostcellsfine_tot_y=nbghostcells+1 118 119 irafx = Agrif_irhox() 120 121 CALL agrif_nemo_init ! specific namelist part if needed 122 123 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),glamt_id) 124 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),glamu_id) 125 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),glamv_id) 126 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),glamf_id) 127 128 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),gphit_id) 129 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),gphiu_id) 130 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),gphiv_id) 131 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),gphif_id) 132 133 ! Horizontal scale factors 134 135 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e1t_id) 136 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e1u_id) 137 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e1v_id) 138 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e1f_id) 139 140 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e2t_id) 141 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e2u_id) 142 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e2v_id) 143 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e2f_id) 144 145 ! Bathymetry 146 147 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),bathy_id) 148 149 ! Vertical scale factors 150 CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3t_id) 151 CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3t_copy_id) 152 CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk+1/),e3t_connect_id) 153 154 CALL agrif_declare_variable((/1,2,0/),(/ind2-1,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3u_id) 155 CALL agrif_declare_variable((/2,1,0/),(/ind2,ind3-1,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3v_id) 156 157 ! Bottom level 158 159 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),bottom_level_id) 160 161 CALL Agrif_Set_bcinterp(glamt_id,interp=AGRIF_linear) 162 CALL Agrif_Set_interp(glamt_id,interp=AGRIF_linear) 163 CALL Agrif_Set_bc( glamt_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 164 165 CALL Agrif_Set_bcinterp(glamu_id,interp=AGRIF_linear) 166 CALL Agrif_Set_interp(glamu_id,interp=AGRIF_linear) 167 CALL Agrif_Set_bc( glamu_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 168 169 CALL Agrif_Set_bcinterp(glamv_id,interp=AGRIF_linear) 170 CALL Agrif_Set_interp(glamv_id,interp=AGRIF_linear) 171 CALL Agrif_Set_bc( glamv_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 172 173 CALL Agrif_Set_bcinterp(glamf_id,interp=AGRIF_linear) 174 CALL Agrif_Set_interp(glamf_id,interp=AGRIF_linear) 175 CALL Agrif_Set_bc( glamf_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 176 177 CALL Agrif_Set_bcinterp(gphit_id,interp=AGRIF_linear) 178 CALL Agrif_Set_interp(gphit_id,interp=AGRIF_linear) 179 CALL Agrif_Set_bc( gphit_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 180 181 CALL Agrif_Set_bcinterp(gphiu_id,interp=AGRIF_linear) 182 CALL Agrif_Set_interp(gphiu_id,interp=AGRIF_linear) 183 CALL Agrif_Set_bc( gphiu_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 184 185 CALL Agrif_Set_bcinterp(gphiv_id,interp=AGRIF_linear) 186 CALL Agrif_Set_interp(gphiv_id,interp=AGRIF_linear) 187 CALL Agrif_Set_bc( gphiv_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 188 189 CALL Agrif_Set_bcinterp(gphif_id,interp=AGRIF_linear) 190 CALL Agrif_Set_interp(gphif_id,interp=AGRIF_linear) 191 CALL Agrif_Set_bc( gphif_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 192 193 ! 194 195 CALL Agrif_Set_bcinterp(e1t_id,interp=AGRIF_ppm) 196 CALL Agrif_Set_interp(e1t_id,interp=AGRIF_ppm) 197 CALL Agrif_Set_bc( e1t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 198 199 CALL Agrif_Set_bcinterp(e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 200 CALL Agrif_Set_interp(e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 201 CALL Agrif_Set_bc( e1u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 202 203 CALL Agrif_Set_bcinterp(e1v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 204 CALL Agrif_Set_interp(e1v_id, interp1=AGRIF_ppm, interp2=Agrif_linear) 205 CALL Agrif_Set_bc( e1v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 206 207 CALL Agrif_Set_bcinterp(e1f_id,interp=AGRIF_linear) 208 CALL Agrif_Set_interp(e1f_id,interp=AGRIF_linear) 209 CALL Agrif_Set_bc( e1f_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 210 211 CALL Agrif_Set_bcinterp(e2t_id,interp=AGRIF_ppm) 212 CALL Agrif_Set_interp(e2t_id,interp=AGRIF_ppm) 213 CALL Agrif_Set_bc( e2t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 214 215 CALL Agrif_Set_bcinterp(e2u_id,interp1=Agrif_linear, interp2=AGRIF_ppm) 216 CALL Agrif_Set_interp(e2u_id,interp1=Agrif_linear, interp2=AGRIF_ppm) 217 CALL Agrif_Set_bc( e2u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 218 219 CALL Agrif_Set_bcinterp(e2v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 220 CALL Agrif_Set_interp(e2v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 221 CALL Agrif_Set_bc( e2v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 222 223 CALL Agrif_Set_bcinterp(e2f_id,interp=AGRIF_linear) 224 CALL Agrif_Set_interp(e2f_id,interp=AGRIF_linear) 225 CALL Agrif_Set_bc( e2f_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 226 227 CALL Agrif_Set_bcinterp(bathy_id,interp=AGRIF_linear) 228 CALL Agrif_Set_interp(bathy_id,interp=AGRIF_linear) 229 CALL Agrif_Set_bc( bathy_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 230 231 ! Vertical scale factors 232 CALL Agrif_Set_bcinterp(e3t_id,interp=AGRIF_ppm) 233 CALL Agrif_Set_interp(e3t_id,interp=AGRIF_ppm) 234 CALL Agrif_Set_bc( e3t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 235 CALL Agrif_Set_Updatetype( e3t_id, update = AGRIF_Update_Average) 236 237 CALL Agrif_Set_bcinterp(e3t_copy_id,interp=AGRIF_constant) 238 CALL Agrif_Set_interp(e3t_copy_id,interp=AGRIF_constant) 239 CALL Agrif_Set_bc( e3t_copy_id, (/-npt_copy*irafx-1,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/)) 240 241 CALL Agrif_Set_bcinterp(e3t_connect_id,interp=AGRIF_ppm) 242 CALL Agrif_Set_interp(e3t_connect_id,interp=AGRIF_ppm) 243 CALL Agrif_Set_bc( e3t_connect_id, (/-(npt_copy+npt_connect)*irafx-1,-npt_copy*irafx-2/)) 244 245 CALL Agrif_Set_bcinterp(e3u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 246 CALL Agrif_Set_interp(e3u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 247 CALL Agrif_Set_bc( e3u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 248 CALL Agrif_Set_Updatetype(e3u_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 249 250 CALL Agrif_Set_bcinterp(e3v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 251 CALL Agrif_Set_interp(e3v_id, interp1=AGRIF_ppm, interp2=Agrif_linear) 252 CALL Agrif_Set_bc( e3v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 253 CALL Agrif_Set_Updatetype(e3v_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 254 255 ! Bottom level 256 CALL Agrif_Set_bcinterp(bottom_level_id,interp=AGRIF_constant) 257 CALL Agrif_Set_interp(bottom_level_id,interp=AGRIF_constant) 258 CALL Agrif_Set_bc( bottom_level_id, (/-npt_copy*irafx-1,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/)) 259 CALL Agrif_Set_Updatetype( bottom_level_id, update = AGRIF_Update_Max) 260 261 end subroutine agrif_declare_var 262 263 264 subroutine agrif_init_lonlat() 265 use agrif_profiles 266 use agrif_util 267 external :: init_glamt, init_glamu, init_glamv, init_glamf 268 external :: init_gphit, init_gphiu, init_gphiv, init_gphif 269 270 call Agrif_Init_variable(glamt_id, procname = init_glamt) 271 call Agrif_Init_variable(glamu_id, procname = init_glamu) 272 call Agrif_Init_variable(glamv_id, procname = init_glamv) 273 call Agrif_Init_variable(glamf_id, procname = init_glamf) 274 275 call Agrif_Init_variable(gphit_id, procname = init_gphit) 276 call Agrif_Init_variable(gphiu_id, procname = init_gphiu) 277 call Agrif_Init_variable(gphiv_id, procname = init_gphiv) 278 call Agrif_Init_variable(gphif_id, procname = init_gphif) 279 280 end subroutine agrif_init_lonlat 281 282 subroutine agrif_init_scales() 283 use agrif_profiles 284 use agrif_util 285 external :: init_e1t, init_e1u, init_e1v, init_e1f 286 external :: init_e2t, init_e2u, init_e2v, init_e2f 287 288 call Agrif_Init_variable(e1t_id, procname = init_e1t) 289 call Agrif_Init_variable(e1u_id, procname = init_e1u) 290 call Agrif_Init_variable(e1v_id, procname = init_e1v) 291 call Agrif_Init_variable(e1f_id, procname = init_e1f) 292 293 call Agrif_Init_variable(e2t_id, procname = init_e2t) 294 call Agrif_Init_variable(e2u_id, procname = init_e2u) 295 call Agrif_Init_variable(e2v_id, procname = init_e2v) 296 call Agrif_Init_variable(e2f_id, procname = init_e2f) 297 298 end subroutine agrif_init_scales 299 300 301 302 SUBROUTINE init_glamt( ptab, i1, i2, j1, j2, before, nb,ndir) 303 use dom_oce 65 66 WRITE(*,*) ' ' 67 WRITE(*,*)'Size of the High resolution grid: ',nx,' x ',ny 68 WRITE(*,*) ' ' 69 70 CALL agrif_init_lonlat() 71 ln_perio = .FALSE. 72 IF( jperio == 1 .OR. jperio == 2 .OR. jperio == 4 ) ln_perio=.TRUE. 73 74 WHERE (glamt < -180) glamt = glamt +360. 75 WHERE (glamt > +180) glamt = glamt -360. 76 77 CALL lbc_lnk( 'glamt', glamt, 'T', 1._wp) 78 CALL lbc_lnk( 'gphit', gphit, 'T', 1._wp) 79 80 WHERE (glamu < -180) glamu = glamu +360. 81 WHERE (glamu > +180) glamu = glamu -360. 82 CALL lbc_lnk( 'glamu', glamu, 'U', 1._wp) 83 CALL lbc_lnk( 'gphiu', gphiu, 'U', 1._wp) 84 85 WHERE (glamv < -180) glamv = glamv +360. 86 WHERE (glamv > +180) glamv = glamv -360. 87 CALL lbc_lnk( 'glamv', glamv, 'V', 1._wp) 88 CALL lbc_lnk( 'gphiv', gphiv, 'V', 1._wp) 89 90 WHERE (glamf < -180) glamf = glamf +360. 91 WHERE (glamf > +180) glamf = glamf -360. 92 CALL lbc_lnk( 'glamf', glamf, 'F', 1._wp) 93 CALL lbc_lnk( 'gphif', gphif, 'F', 1._wp) 94 95 ! Correct South and North 96 IF ((nbondj == -1).OR.(nbondj == 2)) THEN 97 glamt(:,1) = glamt(:,2) 98 gphit(:,1) = gphit(:,2) 99 glamu(:,1) = glamu(:,2) 100 gphiu(:,1) = gphiu(:,2) 101 glamv(:,1) = glamv(:,2) 102 gphiv(:,1) = gphiv(:,2) 103 ENDIF 104 IF ((nbondj == 1).OR.(nbondj == 2)) THEN 105 glamt(:,jpj) = glamt(:,jpj-1) 106 gphit(:,jpj) = gphit(:,jpj-1) 107 glamu(:,jpj) = glamu(:,jpj-1) 108 gphiu(:,jpj) = gphiu(:,jpj-1) 109 glamv(:,jpj) = glamv(:,jpj-1) 110 gphiv(:,jpj) = gphiv(:,jpj-1) 111 glamf(:,jpj) = glamf(:,jpj-1) 112 gphif(:,jpj) = gphif(:,jpj-1) 113 ENDIF 114 115 ! Correct West and East 116 IF( jperio /= 1 ) THEN 117 IF((nbondi == -1) .OR. (nbondi == 2) ) THEN 118 glamt(1,:) = glamt(2,:) 119 gphit(1,:) = gphit(2,:) 120 glamu(1,:) = glamu(2,:) 121 gphiu(1,:) = gphiu(2,:) 122 glamv(1,:) = glamv(2,:) 123 gphiv(1,:) = gphiv(2,:) 124 ENDIF 125 IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 126 glamt(jpi,:) = glamt(jpi-1,:) 127 gphit(jpi,:) = gphit(jpi-1,:) 128 glamu(jpi,:) = glamu(jpi-1,:) 129 gphiu(jpi,:) = gphiu(jpi-1,:) 130 glamv(jpi,:) = glamv(jpi-1,:) 131 gphiv(jpi,:) = gphiv(jpi-1,:) 132 glamf(jpi,:) = glamf(jpi-1,:) 133 gphif(jpi,:) = gphif(jpi-1,:) 134 ENDIF 135 ENDIF 136 137 CALL agrif_init_scales() 138 139 ! Correct South and North 140 IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 141 e1t(:,1) = e1t(:,2) 142 e2t(:,1) = e2t(:,2) 143 e1u(:,1) = e1u(:,2) 144 e2u(:,1) = e2u(:,2) 145 e1v(:,1) = e1v(:,2) 146 e2v(:,1) = e2v(:,2) 147 ENDIF 148 IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 149 e1t(:,jpj) = e1t(:,jpj-1) 150 e2t(:,jpj) = e2t(:,jpj-1) 151 e1u(:,jpj) = e1u(:,jpj-1) 152 e2u(:,jpj) = e2u(:,jpj-1) 153 e1v(:,jpj) = e1v(:,jpj-1) 154 e2v(:,jpj) = e2v(:,jpj-1) 155 e1f(:,jpj) = e1f(:,jpj-1) 156 e2f(:,jpj) = e2f(:,jpj-1) 157 ENDIF 158 159 ! Correct West and East 160 IF( jperio /= 1 ) THEN 161 IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 162 e1t(1,:) = e1t(2,:) 163 e2t(1,:) = e2t(2,:) 164 e1u(1,:) = e1u(2,:) 165 e2u(1,:) = e2u(2,:) 166 e1v(1,:) = e1v(2,:) 167 e2v(1,:) = e2v(2,:) 168 ENDIF 169 IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 170 e1t(jpi,:) = e1t(jpi-1,:) 171 e2t(jpi,:) = e2t(jpi-1,:) 172 e1u(jpi,:) = e1u(jpi-1,:) 173 e2u(jpi,:) = e2u(jpi-1,:) 174 e1v(jpi,:) = e1v(jpi-1,:) 175 e2v(jpi,:) = e2v(jpi-1,:) 176 e1f(jpi,:) = e1f(jpi-1,:) 177 e2f(jpi,:) = e2f(jpi-1,:) 178 ENDIF 179 ENDIF 180 181 END SUBROUTINE Agrif_InitValues_cont 182 183 184 SUBROUTINE agrif_declare_var() 185 !!---------------------------------------------------------------------- 186 !! *** ROUTINE Agrif_InitValues_cont *** 187 !! 188 !! ** Purpose :: Declaration of variables to be interpolated 189 !!---------------------------------------------------------------------- 190 USE par_oce 191 USE dom_oce 192 USE agrif_profiles 193 USE agrif_parameters 194 195 IMPLICIT NONE 196 197 INTEGER :: ind1, ind2, ind3 198 INTEGER :: nx, ny 199 INTEGER ::nbghostcellsfine_tot_x, nbghostcellsfine_tot_y 200 INTEGER :: irafx 201 202 EXTERNAL :: nemo_mapping 203 204 ! 1. Declaration of the type of variable which have to be interpolated 205 !--------------------------------------------------------------------- 206 207 nx=nlci ; ny=nlcj 208 209 ind2 = 2 + nbghostcells_x 210 ind3 = 2 + nbghostcells_y_s 211 nbghostcellsfine_tot_x=nbghostcells_x+1 212 nbghostcellsfine_tot_y=max(nbghostcells_y_s,nbghostcells_y_n)+1 213 214 irafx = Agrif_irhox() 215 216 ! In case of East-West periodicity, prevent AGRIF interpolation at east and west boundaries 217 ! The procnames will not be CALLed at these boundaries 218 if (jperio == 1) THEN 219 CALL Agrif_Set_NearCommonBorderX(.TRUE.) 220 CALL Agrif_Set_DistantCommonBorderX(.TRUE.) 221 endif 222 if (.not.ln_bry_south) THEN 223 CALL Agrif_Set_NearCommonBorderY(.TRUE.) 224 endif 225 226 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),glamt_id) 227 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),glamu_id) 228 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),glamv_id) 229 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),glamf_id) 230 231 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),gphit_id) 232 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),gphiu_id) 233 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),gphiv_id) 234 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),gphif_id) 235 236 ! Horizontal scale factors 237 238 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e1t_id) 239 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e1u_id) 240 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e1v_id) 241 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e1f_id) 242 243 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e2t_id) 244 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e2u_id) 245 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e2v_id) 246 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e2f_id) 247 248 ! Bathymetry 249 250 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),bathy_id) 251 252 ! Vertical scale factors 253 CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3t_id) 254 CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3t_copy_id) 255 CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk+1/),e3t_connect_id) 256 257 CALL agrif_declare_variable((/1,2,0/),(/ind2-1,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3u_id) 258 CALL agrif_declare_variable((/2,1,0/),(/ind2,ind3-1,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3v_id) 259 260 ! Bottom level 261 262 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),bottom_level_id) 263 264 CALL Agrif_Set_bcinterp(glamt_id,interp=AGRIF_linear) 265 CALL Agrif_Set_interp(glamt_id,interp=AGRIF_linear) 266 CALL Agrif_Set_bc( glamt_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 267 268 CALL Agrif_Set_bcinterp(glamu_id,interp=AGRIF_linear) 269 CALL Agrif_Set_interp(glamu_id,interp=AGRIF_linear) 270 CALL Agrif_Set_bc( glamu_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 271 272 CALL Agrif_Set_bcinterp(glamv_id,interp=AGRIF_linear) 273 CALL Agrif_Set_interp(glamv_id,interp=AGRIF_linear) 274 CALL Agrif_Set_bc( glamv_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 275 276 CALL Agrif_Set_bcinterp(glamf_id,interp=AGRIF_linear) 277 CALL Agrif_Set_interp(glamf_id,interp=AGRIF_linear) 278 CALL Agrif_Set_bc( glamf_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 279 280 CALL Agrif_Set_bcinterp(gphit_id,interp=AGRIF_linear) 281 CALL Agrif_Set_interp(gphit_id,interp=AGRIF_linear) 282 CALL Agrif_Set_bc( gphit_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 283 284 CALL Agrif_Set_bcinterp(gphiu_id,interp=AGRIF_linear) 285 CALL Agrif_Set_interp(gphiu_id,interp=AGRIF_linear) 286 CALL Agrif_Set_bc( gphiu_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 287 288 CALL Agrif_Set_bcinterp(gphiv_id,interp=AGRIF_linear) 289 CALL Agrif_Set_interp(gphiv_id,interp=AGRIF_linear) 290 CALL Agrif_Set_bc( gphiv_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 291 292 CALL Agrif_Set_bcinterp(gphif_id,interp=AGRIF_linear) 293 CALL Agrif_Set_interp(gphif_id,interp=AGRIF_linear) 294 CALL Agrif_Set_bc( gphif_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 295 296 ! 297 298 CALL Agrif_Set_bcinterp(e1t_id,interp=AGRIF_ppm) 299 CALL Agrif_Set_interp(e1t_id,interp=AGRIF_ppm) 300 CALL Agrif_Set_bc( e1t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 301 302 CALL Agrif_Set_bcinterp(e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 303 CALL Agrif_Set_interp(e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 304 CALL Agrif_Set_bc( e1u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 305 306 CALL Agrif_Set_bcinterp(e1v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 307 CALL Agrif_Set_interp(e1v_id, interp1=AGRIF_ppm, interp2=Agrif_linear) 308 CALL Agrif_Set_bc( e1v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 309 310 CALL Agrif_Set_bcinterp(e1f_id,interp=AGRIF_linear) 311 CALL Agrif_Set_interp(e1f_id,interp=AGRIF_linear) 312 CALL Agrif_Set_bc( e1f_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 313 314 CALL Agrif_Set_bcinterp(e2t_id,interp=AGRIF_ppm) 315 CALL Agrif_Set_interp(e2t_id,interp=AGRIF_ppm) 316 CALL Agrif_Set_bc( e2t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 317 318 CALL Agrif_Set_bcinterp(e2u_id,interp1=Agrif_linear, interp2=AGRIF_ppm) 319 CALL Agrif_Set_interp(e2u_id,interp1=Agrif_linear, interp2=AGRIF_ppm) 320 CALL Agrif_Set_bc( e2u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 321 322 CALL Agrif_Set_bcinterp(e2v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 323 CALL Agrif_Set_interp(e2v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 324 CALL Agrif_Set_bc( e2v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 325 326 CALL Agrif_Set_bcinterp(e2f_id,interp=AGRIF_linear) 327 CALL Agrif_Set_interp(e2f_id,interp=AGRIF_linear) 328 CALL Agrif_Set_bc( e2f_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 329 330 CALL Agrif_Set_bcinterp(bathy_id,interp=AGRIF_linear) 331 CALL Agrif_Set_interp(bathy_id,interp=AGRIF_linear) 332 CALL Agrif_Set_bc( bathy_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 333 334 ! Vertical scale factors 335 CALL Agrif_Set_bcinterp(e3t_id,interp=AGRIF_ppm) 336 CALL Agrif_Set_interp(e3t_id,interp=AGRIF_ppm) 337 CALL Agrif_Set_bc( e3t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 338 CALL Agrif_Set_Updatetype( e3t_id, update = AGRIF_Update_Average) 339 340 CALL Agrif_Set_bcinterp(e3t_copy_id,interp=AGRIF_constant) 341 CALL Agrif_Set_interp(e3t_copy_id,interp=AGRIF_constant) 342 CALL Agrif_Set_bc( e3t_copy_id, (/-npt_copy*irafx,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/)) 343 344 CALL Agrif_Set_bcinterp(e3t_connect_id,interp=AGRIF_ppm) 345 CALL Agrif_Set_interp(e3t_connect_id,interp=AGRIF_ppm) 346 CALL Agrif_Set_bc( e3t_connect_id, (/-(npt_copy+npt_connect)*irafx-1,-npt_copy*irafx/)) 347 348 CALL Agrif_Set_bcinterp(e3u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 349 CALL Agrif_Set_interp(e3u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 350 CALL Agrif_Set_bc( e3u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 351 CALL Agrif_Set_Updatetype(e3u_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 352 353 CALL Agrif_Set_bcinterp(e3v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 354 CALL Agrif_Set_interp(e3v_id, interp1=AGRIF_ppm, interp2=Agrif_linear) 355 CALL Agrif_Set_bc( e3v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 356 CALL Agrif_Set_Updatetype(e3v_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 357 358 ! Bottom level 359 CALL Agrif_Set_bcinterp(bottom_level_id,interp=AGRIF_constant) 360 CALL Agrif_Set_interp(bottom_level_id,interp=AGRIF_constant) 361 CALL Agrif_Set_bc( bottom_level_id, (/-npt_copy*irafx,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/)) 362 CALL Agrif_Set_Updatetype( bottom_level_id, update = AGRIF_Update_Max) 363 364 CALL Agrif_Set_ExternalMapping(nemo_mapping) 365 366 END SUBROUTINE agrif_declare_var 367 368 SUBROUTINE nemo_mapping(ndim,ptx,pty,bounds,bounds_chunks,correction_required,nb_chunks) 369 USE dom_oce 370 INTEGER :: ndim 371 INTEGER :: ptx, pty 372 INTEGER,DIMENSION(ndim,2,2) :: bounds 373 INTEGER,DIMENSION(:,:,:,:),allocatable :: bounds_chunks 374 LOGICAL,DIMENSION(:),allocatable :: correction_required 375 INTEGER :: nb_chunks 376 INTEGER :: i 377 378 IF (agrif_debug_interp) THEN 379 DO i = 1, ndim 380 print *,'direction = ',i,bounds(i,1,2),bounds(i,2,2) 381 END DO 382 ENDIF 383 384 IF( bounds(2,2,2) > jpjglo ) THEN 385 IF( bounds(2,1,2) <= jpjglo ) THEN 386 nb_chunks = 2 387 ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2)) 388 ALLOCATE(correction_required(nb_chunks)) 389 DO i = 1, nb_chunks 390 bounds_chunks(i,:,:,:) = bounds 391 END DO 392 393 ! FIRST CHUNCK (for j<=jpjglo) 394 395 ! Original indices 396 bounds_chunks(1,1,1,1) = bounds(1,1,2) 397 bounds_chunks(1,1,2,1) = bounds(1,2,2) 398 bounds_chunks(1,2,1,1) = bounds(2,1,2) 399 bounds_chunks(1,2,2,1) = jpjglo 400 401 bounds_chunks(1,1,1,2) = bounds(1,1,2) 402 bounds_chunks(1,1,2,2) = bounds(1,2,2) 403 bounds_chunks(1,2,1,2) = bounds(2,1,2) 404 bounds_chunks(1,2,2,2) = jpjglo 405 406 ! Correction required or not 407 correction_required(1)=.FALSE. 408 409 ! SECOND CHUNCK (for j>jpjglo) 410 411 !Original indices 412 bounds_chunks(2,1,1,1) = bounds(1,1,2) 413 bounds_chunks(2,1,2,1) = bounds(1,2,2) 414 bounds_chunks(2,2,1,1) = jpjglo-2 415 bounds_chunks(2,2,2,1) = bounds(2,2,2) 416 417 ! Where to find them 418 ! We use the relation TAB(ji,jj)=TAB(jpiglo-ji+2,jpjglo-2-(jj-jpjglo)) 419 420 IF (ptx == 2) THEN ! T, V points 421 bounds_chunks(2,1,1,2) = jpiglo-bounds(1,2,2)+2 422 bounds_chunks(2,1,2,2) = jpiglo-bounds(1,1,2)+2 423 ELSE ! U, F points 424 bounds_chunks(2,1,1,2) = jpiglo-bounds(1,2,2)+1 425 bounds_chunks(2,1,2,2) = jpiglo-bounds(1,1,2)+1 426 ENDIF 427 428 IF (pty == 2) THEN ! T, U points 429 bounds_chunks(2,2,1,2) = jpjglo-2-(bounds(2,2,2) -jpjglo) 430 bounds_chunks(2,2,2,2) = jpjglo-2-(jpjglo-2 -jpjglo) 431 ELSE ! V, F points 432 bounds_chunks(2,2,1,2) = jpjglo-3-(bounds(2,2,2) -jpjglo) 433 bounds_chunks(2,2,2,2) = jpjglo-3-(jpjglo-2 -jpjglo) 434 ENDIF 435 436 ! Correction required or not 437 correction_required(2)=.TRUE. 438 439 ELSE 440 441 nb_chunks = 1 442 ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2)) 443 ALLOCATE(correction_required(nb_chunks)) 444 DO i=1,nb_chunks 445 bounds_chunks(i,:,:,:) = bounds 446 END DO 447 448 bounds_chunks(1,1,1,1) = bounds(1,1,2) 449 bounds_chunks(1,1,2,1) = bounds(1,2,2) 450 bounds_chunks(1,2,1,1) = bounds(2,1,2) 451 bounds_chunks(1,2,2,1) = bounds(2,2,2) 452 453 bounds_chunks(1,1,1,2) = jpiglo-bounds(1,2,2)+2 454 bounds_chunks(1,1,2,2) = jpiglo-bounds(1,1,2)+2 455 456 bounds_chunks(1,2,1,2) = jpjglo-2-(bounds(2,2,2)-jpjglo) 457 bounds_chunks(1,2,2,2) = jpjglo-2-(bounds(2,1,2)-jpjglo) 458 459 IF (ptx == 2) THEN ! T, V points 460 bounds_chunks(1,1,1,2) = jpiglo-bounds(1,2,2)+2 461 bounds_chunks(1,1,2,2) = jpiglo-bounds(1,1,2)+2 462 ELSE ! U, F points 463 bounds_chunks(1,1,1,2) = jpiglo-bounds(1,2,2)+1 464 bounds_chunks(1,1,2,2) = jpiglo-bounds(1,1,2)+1 465 ENDIF 466 467 IF (pty == 2) THEN ! T, U points 468 bounds_chunks(1,2,1,2) = jpjglo-2-(bounds(2,2,2) -jpjglo) 469 bounds_chunks(1,2,2,2) = jpjglo-2-(bounds(2,1,2) -jpjglo) 470 ELSE ! V, F points 471 bounds_chunks(1,2,1,2) = jpjglo-3-(bounds(2,2,2) -jpjglo) 472 bounds_chunks(1,2,2,2) = jpjglo-3-(bounds(2,1,2) -jpjglo) 473 ENDIF 474 475 correction_required(1)=.TRUE. 476 477 ENDIF ! bounds(2,1,2) <= jpjglo 478 479 ELSE IF (bounds(1,1,2) < 1) THEN 480 481 IF (bounds(1,2,2) > 0) THEN 482 nb_chunks = 2 483 ALLOCATE(correction_required(nb_chunks)) 484 correction_required=.FALSE. 485 ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2)) 486 DO i=1,nb_chunks 487 bounds_chunks(i,:,:,:) = bounds 488 END DO 489 490 bounds_chunks(1,1,1,2) = bounds(1,1,2)+jpiglo-2 491 bounds_chunks(1,1,2,2) = 1+jpiglo-2 492 493 bounds_chunks(1,1,1,1) = bounds(1,1,2) 494 bounds_chunks(1,1,2,1) = 1 495 496 bounds_chunks(2,1,1,2) = 2 497 bounds_chunks(2,1,2,2) = bounds(1,2,2) 498 499 bounds_chunks(2,1,1,1) = 2 500 bounds_chunks(2,1,2,1) = bounds(1,2,2) 501 ELSE 502 nb_chunks = 1 503 ALLOCATE(correction_required(nb_chunks)) 504 correction_required=.FALSE. 505 ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2)) 506 DO i=1,nb_chunks 507 bounds_chunks(i,:,:,:) = bounds 508 END DO 509 bounds_chunks(1,1,1,2) = bounds(1,1,2)+jpiglo-2 510 bounds_chunks(1,1,2,2) = bounds(1,2,2)+jpiglo-2 511 512 bounds_chunks(1,1,1,1) = bounds(1,1,2) 513 bounds_chunks(1,1,2,1) = bounds(1,2,2) 514 ENDIF 515 516 ELSE 517 518 nb_chunks=1 519 ALLOCATE(correction_required(nb_chunks)) 520 correction_required=.FALSE. 521 ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2)) 522 DO i=1,nb_chunks 523 bounds_chunks(i,:,:,:) = bounds 524 END DO 525 bounds_chunks(1,1,1,2) = bounds(1,1,2) 526 bounds_chunks(1,1,2,2) = bounds(1,2,2) 527 bounds_chunks(1,2,1,2) = bounds(2,1,2) 528 bounds_chunks(1,2,2,2) = bounds(2,2,2) 529 530 bounds_chunks(1,1,1,1) = bounds(1,1,2) 531 bounds_chunks(1,1,2,1) = bounds(1,2,2) 532 bounds_chunks(1,2,1,1) = bounds(2,1,2) 533 bounds_chunks(1,2,2,1) = bounds(2,2,2) 534 535 ENDIF 536 537 END SUBROUTINE nemo_mapping 538 539 FUNCTION agrif_external_switch_index(ptx,pty,i1,isens) 540 USE dom_oce 541 INTEGER :: ptx, pty, i1, isens 542 INTEGER :: agrif_external_switch_index 543 544 IF( isens == 1 ) THEN 545 IF( ptx == 2 ) THEN ! T, V points 546 agrif_external_switch_index = jpiglo-i1+2 547 ELSE ! U, F points 548 agrif_external_switch_index = jpiglo-i1+1 549 ENDIF 550 ELSE IF (isens ==2) THEN 551 IF (pty == 2) THEN ! T, U points 552 agrif_external_switch_index = jpjglo-2-(i1 -jpjglo) 553 ELSE ! V, F points 554 agrif_external_switch_index = jpjglo-3-(i1 -jpjglo) 555 ENDIF 556 ENDIF 557 558 END FUNCTION agrif_external_switch_index 559 560 SUBROUTINE correct_field(tab2d,i1,i2,j1,j2) 561 USE dom_oce 562 INTEGER :: i1,i2,j1,j2 563 REAL,DIMENSION(i1:i2,j1:j2) :: tab2d 564 565 INTEGER :: i,j 566 REAL,DIMENSION(i1:i2,j1:j2) :: tab2dtemp 567 568 tab2dtemp = tab2d 569 570 DO j=j1,j2 571 DO i=i1,i2 572 tab2d(i,j)=tab2dtemp(i2-(i-i1),j2-(j-j1)) 573 END DO 574 ENDDO 575 576 END SUBROUTINE correct_field 577 578 SUBROUTINE agrif_init_lonlat() 579 USE agrif_profiles 580 USE agrif_util 581 USE dom_oce 582 583 EXTERNAL :: init_glamt, init_glamu, init_glamv, init_glamf 584 EXTERNAL :: init_gphit, init_gphiu, init_gphiv, init_gphif 585 EXTERNAL :: longitude_linear_interp 586 587 INTEGER :: ji,jj,i1,i2,j1,j2 588 REAL, DIMENSION(jpi,jpj) :: tab2dtemp 589 INTEGER :: ind2,ind3 590 INTEGER :: irhox, irhoy 591 592 irhox = agrif_irhox() 593 irhoy = agrif_irhoy() 594 CALL Agrif_Set_external_linear_interp(longitude_linear_interp) 595 596 CALL Agrif_Init_variable(glamt_id, procname = init_glamt) 597 CALL Agrif_Init_variable(glamu_id, procname = init_glamu) 598 CALL Agrif_Init_variable(glamv_id, procname = init_glamv) 599 CALL Agrif_Init_variable(glamf_id, procname = init_glamf) 600 CALL Agrif_UnSet_external_linear_interp() 601 602 CALL Agrif_Init_variable(gphit_id, procname = init_gphit) 603 CALL Agrif_Init_variable(gphiu_id, procname = init_gphiu) 604 CALL Agrif_Init_variable(gphiv_id, procname = init_gphiv) 605 CALL Agrif_Init_variable(gphif_id, procname = init_gphif) 606 607 END SUBROUTINE agrif_init_lonlat 608 609 REAL FUNCTION longitude_linear_interp(x1,x2,coeff) 610 REAL :: x1, x2, coeff 611 REAL :: val_interp 612 613 IF( (x1*x2 <= -50*50) ) THEN 614 IF( x1 < 0 ) THEN 615 val_interp = coeff *(x1+360.) + (1.-coeff) *x2 616 ELSE 617 val_interp = coeff *x1 + (1.-coeff) *(x2+360.) 618 ENDIF 619 IF ((val_interp) >=180.) val_interp = val_interp - 360. 620 ELSE 621 val_interp = coeff * x1 + (1.-coeff) * x2 622 ENDIF 623 longitude_linear_interp = val_interp 624 625 END FUNCTION longitude_linear_interp 626 627 SUBROUTINE agrif_init_scales() 628 USE agrif_profiles 629 USE agrif_util 630 USE dom_oce 631 USE lbclnk 632 LOGICAL :: ln_perio 633 INTEGER nx,ny 634 635 EXTERNAL :: init_e1t, init_e1u, init_e1v, init_e1f 636 EXTERNAL :: init_e2t, init_e2u, init_e2v, init_e2f 637 638 nx = nlci ; ny = nlcj 639 640 ln_perio=.FALSE. 641 if( jperio ==1 .OR. jperio==2 .OR. jperio==4) ln_perio=.TRUE. 642 643 CALL Agrif_Init_variable(e1t_id, procname = init_e1t) 644 CALL Agrif_Init_variable(e1u_id, procname = init_e1u) 645 CALL Agrif_Init_variable(e1v_id, procname = init_e1v) 646 CALL Agrif_Init_variable(e1f_id, procname = init_e1f) 647 648 CALL Agrif_Init_variable(e2t_id, procname = init_e2t) 649 CALL Agrif_Init_variable(e2u_id, procname = init_e2u) 650 CALL Agrif_Init_variable(e2v_id, procname = init_e2v) 651 CALL Agrif_Init_variable(e2f_id, procname = init_e2f) 652 653 CALL lbc_lnk( 'e1t', e1t, 'T', 1._wp) 654 CALL lbc_lnk( 'e2t', e2t, 'T', 1._wp) 655 CALL lbc_lnk( 'e1u', e1u, 'U', 1._wp) 656 CALL lbc_lnk( 'e2u', e2u, 'U', 1._wp) 657 CALL lbc_lnk( 'e1v', e1v, 'V', 1._wp) 658 CALL lbc_lnk( 'e2v', e2v, 'V', 1._wp) 659 CALL lbc_lnk( 'e1f', e1f, 'F', 1._wp) 660 CALL lbc_lnk( 'e2f', e2f, 'F', 1._wp) 661 662 END SUBROUTINE agrif_init_scales 663 664 SUBROUTINE init_glamt( ptab, i1, i2, j1, j2, before, nb,ndir) 665 USE dom_oce 304 666 !!---------------------------------------------------------------------- 305 667 !! *** ROUTINE interpsshn *** 306 !!---------------------------------------------------------------------- 668 !!---------------------------------------------------------------------- 669 INTEGER , INTENT(in ) :: i1, i2, j1, j2 670 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 671 LOGICAL , INTENT(in ) :: before 672 INTEGER , INTENT(in ) :: nb , ndir 673 674 ! 675 !!---------------------------------------------------------------------- 676 ! 677 IF( before) THEN 678 ptab(i1:i2,j1:j2) = glamt(i1:i2,j1:j2) 679 ELSE 680 glamt(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 681 ENDIF 682 ! 683 END SUBROUTINE init_glamt 684 685 SUBROUTINE init_glamu( ptab, i1, i2, j1, j2, before, nb,ndir) 686 USE dom_oce 687 !!---------------------------------------------------------------------- 688 !! *** ROUTINE interpsshn *** 689 !!---------------------------------------------------------------------- 307 690 INTEGER , INTENT(in ) :: i1, i2, j1, j2 308 691 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab … … 311 694 LOGICAL :: western_side, eastern_side,northern_side,southern_side 312 695 ! 313 !!---------------------------------------------------------------------- 314 ! 315 western_side = (nb == 1).AND.(ndir == 1) 316 eastern_side = (nb == 1).AND.(ndir == 2) 317 southern_side = (nb == 2).AND.(ndir == 1) 318 northern_side = (nb == 2).AND.(ndir == 2) 319 IF( before) THEN 320 ptab(i1:i2,j1:j2) = glamt(i1:i2,j1:j2) 321 ELSE 322 glamt(i1:i2,j1:j2)=ptab 323 ENDIF 324 ! 325 END SUBROUTINE init_glamt 326 327 SUBROUTINE init_glamu( ptab, i1, i2, j1, j2, before, nb,ndir) 328 use dom_oce 696 !!---------------------------------------------------------------------- 697 ! 698 IF( before) THEN 699 ptab(i1:i2,j1:j2) = glamu(i1:i2,j1:j2) 700 ELSE 701 glamu(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 702 ENDIF 703 ! 704 END SUBROUTINE init_glamu 705 706 SUBROUTINE init_glamv( ptab, i1, i2, j1, j2, before, nb,ndir) 707 USE dom_oce 329 708 !!---------------------------------------------------------------------- 330 709 !! *** ROUTINE interpsshn *** 331 !!---------------------------------------------------------------------- 332 INTEGER , INTENT(in ) :: i1, i2, j1, j2 333 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 334 LOGICAL , INTENT(in ) :: before 335 INTEGER , INTENT(in ) :: nb , ndir 336 LOGICAL :: western_side, eastern_side,northern_side,southern_side 337 ! 338 !!---------------------------------------------------------------------- 339 ! 340 western_side = (nb == 1).AND.(ndir == 1) 341 eastern_side = (nb == 1).AND.(ndir == 2) 342 southern_side = (nb == 2).AND.(ndir == 1) 343 northern_side = (nb == 2).AND.(ndir == 2) 344 IF( before) THEN 345 ptab(i1:i2,j1:j2) = glamu(i1:i2,j1:j2) 346 ELSE 347 glamu(i1:i2,j1:j2)=ptab 348 ENDIF 349 ! 350 END SUBROUTINE init_glamu 351 352 SUBROUTINE init_glamv( ptab, i1, i2, j1, j2, before, nb,ndir) 353 use dom_oce 710 !!---------------------------------------------------------------------- 711 INTEGER , INTENT(in ) :: i1, i2, j1, j2 712 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 713 LOGICAL , INTENT(in ) :: before 714 INTEGER , INTENT(in ) :: nb , ndir 715 ! 716 !!---------------------------------------------------------------------- 717 ! 718 IF( before) THEN 719 ptab(i1:i2,j1:j2) = glamv(i1:i2,j1:j2) 720 ELSE 721 glamv(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 722 ENDIF 723 ! 724 END SUBROUTINE init_glamv 725 726 SUBROUTINE init_glamf( ptab, i1, i2, j1, j2, before, nb,ndir) 727 USE dom_oce 728 !!---------------------------------------------------------------------- 729 !! *** ROUTINE init_glamf *** 730 !!---------------------------------------------------------------------- 731 INTEGER , INTENT(in ) :: i1, i2, j1, j2 732 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 733 LOGICAL , INTENT(in ) :: before 734 INTEGER , INTENT(in ) :: nb , ndir 735 ! 736 !!---------------------------------------------------------------------- 737 ! 738 IF( before) THEN 739 ptab(i1:i2,j1:j2) = glamf(i1:i2,j1:j2) 740 ELSE 741 glamf(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 742 ENDIF 743 ! 744 END SUBROUTINE init_glamf 745 746 SUBROUTINE init_gphit( ptab, i1, i2, j1, j2, before, nb,ndir) 747 USE dom_oce 748 !!---------------------------------------------------------------------- 749 !! *** ROUTINE init_gphit *** 750 !!---------------------------------------------------------------------- 751 INTEGER , INTENT(in ) :: i1, i2, j1, j2 752 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 753 LOGICAL , INTENT(in ) :: before 754 INTEGER , INTENT(in ) :: nb , ndir 755 ! 756 !!---------------------------------------------------------------------- 757 ! 758 IF( before) THEN 759 ptab(i1:i2,j1:j2) = gphit(i1:i2,j1:j2) 760 ELSE 761 gphit(i1:i2,j1:j2)=ptab 762 ENDIF 763 ! 764 END SUBROUTINE init_gphit 765 766 SUBROUTINE init_gphiu( ptab, i1, i2, j1, j2, before, nb,ndir) 767 USE dom_oce 768 !!---------------------------------------------------------------------- 769 !! *** ROUTINE init_gphiu *** 770 !!---------------------------------------------------------------------- 771 INTEGER , INTENT(in ) :: i1, i2, j1, j2 772 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 773 LOGICAL , INTENT(in ) :: before 774 INTEGER , INTENT(in ) :: nb , ndir 775 ! 776 !!---------------------------------------------------------------------- 777 ! 778 IF( before) THEN 779 ptab(i1:i2,j1:j2) = gphiu(i1:i2,j1:j2) 780 ELSE 781 gphiu(i1:i2,j1:j2)=ptab 782 ENDIF 783 ! 784 END SUBROUTINE init_gphiu 785 786 SUBROUTINE init_gphiv( ptab, i1, i2, j1, j2, before, nb,ndir) 787 USE dom_oce 788 !!---------------------------------------------------------------------- 789 !! *** ROUTINE init_gphiv *** 790 !!---------------------------------------------------------------------- 791 INTEGER , INTENT(in ) :: i1, i2, j1, j2 792 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 793 LOGICAL , INTENT(in ) :: before 794 INTEGER , INTENT(in ) :: nb , ndir 795 ! 796 !!---------------------------------------------------------------------- 797 798 IF( before) THEN 799 ptab(i1:i2,j1:j2) = gphiv(i1:i2,j1:j2) 800 ELSE 801 gphiv(i1:i2,j1:j2)=ptab 802 ENDIF 803 ! 804 END SUBROUTINE init_gphiv 805 806 807 SUBROUTINE init_gphif( ptab, i1, i2, j1, j2, before, nb,ndir) 808 USE dom_oce 809 !!---------------------------------------------------------------------- 810 !! *** ROUTINE init_gphif *** 811 !!---------------------------------------------------------------------- 812 INTEGER , INTENT(in ) :: i1, i2, j1, j2 813 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 814 LOGICAL , INTENT(in ) :: before 815 INTEGER , INTENT(in ) :: nb , ndir 816 ! 817 !!---------------------------------------------------------------------- 818 ! 819 IF( before) THEN 820 ptab(i1:i2,j1:j2) = gphif(i1:i2,j1:j2) 821 ELSE 822 gphif(i1:i2,j1:j2)=ptab 823 ENDIF 824 ! 825 END SUBROUTINE init_gphif 826 827 828 SUBROUTINE init_e1t( ptab, i1, i2, j1, j2, before, nb,ndir) 829 USE dom_oce 830 USE agrif_parameters 831 !!---------------------------------------------------------------------- 832 !! *** ROUTINE init_e1t *** 833 !!---------------------------------------------------------------------- 834 INTEGER , INTENT(in ) :: i1, i2, j1, j2 835 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 836 LOGICAL , INTENT(in ) :: before 837 INTEGER , INTENT(in ) :: nb , ndir 838 ! 839 !!---------------------------------------------------------------------- 840 ! 841 INTEGER :: jj 842 843 IF( before) THEN 844 ! May need to extend at south boundary 845 IF (j1<1) THEN 846 IF (.NOT.agrif_child(ln_bry_south)) THEN 847 IF ((nbondj == -1).OR.(nbondj == 2)) THEN 848 DO jj=1,j2 849 ptab(i1:i2,jj)=e1t(i1:i2,jj) 850 ENDDO 851 DO jj=j1,0 852 ptab(i1:i2,jj)=e1t(i1:i2,1) 853 ENDDO 854 ENDIF 855 ELSE 856 stop "OUT OF BOUNDS" 857 ENDIF 858 ELSE 859 ptab(i1:i2,j1:j2) = e1t(i1:i2,j1:j2) 860 ENDIF 861 ELSE 862 e1t(i1:i2,j1:j2)=ptab/Agrif_rhoy() 863 ENDIF 864 ! 865 END SUBROUTINE init_e1t 866 867 SUBROUTINE init_e1u( ptab, i1, i2, j1, j2, before, nb,ndir) 868 USE dom_oce 869 USE agrif_parameters 870 !!---------------------------------------------------------------------- 871 !! *** ROUTINE init_e1u *** 872 !!---------------------------------------------------------------------- 873 INTEGER , INTENT(in ) :: i1, i2, j1, j2 874 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 875 LOGICAL , INTENT(in ) :: before 876 INTEGER , INTENT(in ) :: nb , ndir 877 ! 878 !!---------------------------------------------------------------------- 879 ! 880 INTEGER :: jj 881 882 IF( before) THEN 883 IF (j1<1) THEN 884 IF (.NOT.agrif_child(ln_bry_south)) THEN 885 IF ((nbondj == -1).OR.(nbondj == 2)) THEN 886 DO jj=1,j2 887 ptab(i1:i2,jj)=e1u(i1:i2,jj) 888 ENDDO 889 DO jj=j1,0 890 ptab(i1:i2,jj)=e1u(i1:i2,1) 891 ENDDO 892 ENDIF 893 ELSE 894 stop "OUT OF BOUNDS" 895 ENDIF 896 ELSE 897 ptab(i1:i2,j1:j2) = e1u(i1:i2,j1:j2) 898 ENDIF 899 ELSE 900 e1u(i1:i2,j1:j2)=ptab/Agrif_rhoy() 901 ENDIF 902 ! 903 END SUBROUTINE init_e1u 904 905 SUBROUTINE init_e1v( ptab, i1, i2, j1, j2, before, nb,ndir) 906 USE dom_oce 907 !!---------------------------------------------------------------------- 908 !! *** ROUTINE init_e1v *** 909 !!---------------------------------------------------------------------- 910 INTEGER , INTENT(in ) :: i1, i2, j1, j2 911 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 912 LOGICAL , INTENT(in ) :: before 913 INTEGER , INTENT(in ) :: nb , ndir 914 ! 915 !!---------------------------------------------------------------------- 916 ! 917 IF( before) THEN 918 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) 919 ELSE 920 e1v(i1:i2,j1:j2)=ptab/Agrif_rhoy() 921 ENDIF 922 ! 923 END SUBROUTINE init_e1v 924 925 SUBROUTINE init_e1f( ptab, i1, i2, j1, j2, before, nb,ndir) 926 USE dom_oce 927 !!---------------------------------------------------------------------- 928 !! *** ROUTINE init_e1f *** 929 !!---------------------------------------------------------------------- 930 INTEGER , INTENT(in ) :: i1, i2, j1, j2 931 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 932 LOGICAL , INTENT(in ) :: before 933 INTEGER , INTENT(in ) :: nb , ndir 934 ! 935 !!---------------------------------------------------------------------- 936 ! 937 IF( before) THEN 938 ptab(i1:i2,j1:j2) = e1f(i1:i2,j1:j2) 939 ELSE 940 e1f(i1:i2,j1:j2)=ptab/Agrif_rhoy() 941 ENDIF 942 ! 943 END SUBROUTINE init_e1f 944 945 SUBROUTINE init_e2t( ptab, i1, i2, j1, j2, before, nb,ndir) 946 USE dom_oce 947 USE agrif_parameters 948 !!---------------------------------------------------------------------- 949 !! *** ROUTINE init_e2t *** 950 !!---------------------------------------------------------------------- 951 INTEGER , INTENT(in ) :: i1, i2, j1, j2 952 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 953 LOGICAL , INTENT(in ) :: before 954 INTEGER , INTENT(in ) :: nb , ndir 955 ! 956 !!---------------------------------------------------------------------- 957 ! 958 INTEGER :: jj 959 960 IF( before) THEN 961 IF (j1<1) THEN 962 IF (.NOT.agrif_child(ln_bry_south)) THEN 963 IF ((nbondj == -1).OR.(nbondj == 2)) THEN 964 DO jj=1,j2 965 ptab(i1:i2,jj)=e2t(i1:i2,jj) 966 ENDDO 967 DO jj=j1,0 968 ptab(i1:i2,jj)=e2t(i1:i2,1) 969 ENDDO 970 ENDIF 971 ELSE 972 stop "OUT OF BOUNDS" 973 ENDIF 974 ELSE 975 ptab(i1:i2,j1:j2) = e2t(i1:i2,j1:j2) 976 ENDIF 977 ELSE 978 e2t(i1:i2,j1:j2)=ptab/Agrif_rhoy() 979 ENDIF 980 ! 981 END SUBROUTINE init_e2t 982 983 SUBROUTINE init_e2u( ptab, i1, i2, j1, j2, before, nb,ndir) 984 USE dom_oce 985 USE agrif_parameters 354 986 !!---------------------------------------------------------------------- 355 987 !! *** ROUTINE interpsshn *** 356 !!---------------------------------------------------------------------- 357 INTEGER , INTENT(in ) :: i1, i2, j1, j2 358 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 359 LOGICAL , INTENT(in ) :: before 360 INTEGER , INTENT(in ) :: nb , ndir 361 LOGICAL :: western_side, eastern_side,northern_side,southern_side 362 ! 363 !!---------------------------------------------------------------------- 364 ! 365 western_side = (nb == 1).AND.(ndir == 1) 366 eastern_side = (nb == 1).AND.(ndir == 2) 367 southern_side = (nb == 2).AND.(ndir == 1) 368 northern_side = (nb == 2).AND.(ndir == 2) 369 IF( before) THEN 370 ptab(i1:i2,j1:j2) = glamv(i1:i2,j1:j2) 371 ELSE 372 glamv(i1:i2,j1:j2)=ptab 373 ENDIF 374 ! 375 END SUBROUTINE init_glamv 376 377 SUBROUTINE init_glamf( ptab, i1, i2, j1, j2, before, nb,ndir) 378 use dom_oce 988 !!---------------------------------------------------------------------- 989 INTEGER , INTENT(in ) :: i1, i2, j1, j2 990 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 991 LOGICAL , INTENT(in ) :: before 992 INTEGER , INTENT(in ) :: nb , ndir 993 ! 994 !!---------------------------------------------------------------------- 995 ! 996 INTEGER :: jj 997 998 IF( before) THEN 999 IF (j1<1) THEN 1000 IF (.NOT.agrif_child(ln_bry_south)) THEN 1001 IF ((nbondj == -1).OR.(nbondj == 2)) THEN 1002 DO jj=1,j2 1003 ptab(i1:i2,jj)=e2u(i1:i2,jj) 1004 ENDDO 1005 DO jj=j1,0 1006 ptab(i1:i2,jj)=e2u(i1:i2,1) 1007 ENDDO 1008 ENDIF 1009 ELSE 1010 stop "OUT OF BOUNDS" 1011 ENDIF 1012 ELSE 1013 ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) 1014 ENDIF 1015 ELSE 1016 e2u(i1:i2,j1:j2)=ptab/Agrif_rhoy() 1017 ENDIF 1018 ! 1019 END SUBROUTINE init_e2u 1020 1021 SUBROUTINE init_e2v( ptab, i1, i2, j1, j2, before, nb,ndir) 1022 USE dom_oce 379 1023 !!---------------------------------------------------------------------- 380 1024 !! *** ROUTINE interpsshn *** 381 !!---------------------------------------------------------------------- 382 INTEGER , INTENT(in ) :: i1, i2, j1, j2 383 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 384 LOGICAL , INTENT(in ) :: before 385 INTEGER , INTENT(in ) :: nb , ndir 386 LOGICAL :: western_side, eastern_side,northern_side,southern_side 387 ! 388 !!---------------------------------------------------------------------- 389 ! 390 western_side = (nb == 1).AND.(ndir == 1) 391 eastern_side = (nb == 1).AND.(ndir == 2) 392 southern_side = (nb == 2).AND.(ndir == 1) 393 northern_side = (nb == 2).AND.(ndir == 2) 394 IF( before) THEN 395 ptab(i1:i2,j1:j2) = glamf(i1:i2,j1:j2) 396 ELSE 397 glamf(i1:i2,j1:j2)=ptab 398 ENDIF 399 ! 400 END SUBROUTINE init_glamf 401 402 SUBROUTINE init_gphit( ptab, i1, i2, j1, j2, before, nb,ndir) 403 use dom_oce 1025 !!---------------------------------------------------------------------- 1026 INTEGER , INTENT(in ) :: i1, i2, j1, j2 1027 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1028 LOGICAL , INTENT(in ) :: before 1029 INTEGER , INTENT(in ) :: nb , ndir 1030 ! 1031 !!---------------------------------------------------------------------- 1032 IF( before) THEN 1033 ptab(i1:i2,j1:j2) = e2v(i1:i2,j1:j2) 1034 ELSE 1035 e2v(i1:i2,j1:j2)=ptab/Agrif_rhoy() 1036 ENDIF 1037 ! 1038 END SUBROUTINE init_e2v 1039 1040 SUBROUTINE init_e2f( ptab, i1, i2, j1, j2, before, nb,ndir) 1041 USE dom_oce 404 1042 !!---------------------------------------------------------------------- 405 1043 !! *** ROUTINE interpsshn *** 406 !!---------------------------------------------------------------------- 407 INTEGER , INTENT(in ) :: i1, i2, j1, j2 408 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 409 LOGICAL , INTENT(in ) :: before 410 INTEGER , INTENT(in ) :: nb , ndir 411 LOGICAL :: western_side, eastern_side,northern_side,southern_side 412 ! 413 !!---------------------------------------------------------------------- 414 ! 415 western_side = (nb == 1).AND.(ndir == 1) 416 eastern_side = (nb == 1).AND.(ndir == 2) 417 southern_side = (nb == 2).AND.(ndir == 1) 418 northern_side = (nb == 2).AND.(ndir == 2) 419 IF( before) THEN 420 ptab(i1:i2,j1:j2) = gphit(i1:i2,j1:j2) 421 ELSE 422 gphit(i1:i2,j1:j2)=ptab 423 ENDIF 424 ! 425 END SUBROUTINE init_gphit 426 427 SUBROUTINE init_gphiu( ptab, i1, i2, j1, j2, before, nb,ndir) 428 use dom_oce 429 !!---------------------------------------------------------------------- 430 !! *** ROUTINE interpsshn *** 431 !!---------------------------------------------------------------------- 432 INTEGER , INTENT(in ) :: i1, i2, j1, j2 433 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 434 LOGICAL , INTENT(in ) :: before 435 INTEGER , INTENT(in ) :: nb , ndir 436 LOGICAL :: western_side, eastern_side,northern_side,southern_side 437 ! 438 !!---------------------------------------------------------------------- 439 ! 440 western_side = (nb == 1).AND.(ndir == 1) 441 eastern_side = (nb == 1).AND.(ndir == 2) 442 southern_side = (nb == 2).AND.(ndir == 1) 443 northern_side = (nb == 2).AND.(ndir == 2) 444 IF( before) THEN 445 ptab(i1:i2,j1:j2) = gphiu(i1:i2,j1:j2) 446 ELSE 447 gphiu(i1:i2,j1:j2)=ptab 448 ENDIF 449 ! 450 END SUBROUTINE init_gphiu 451 452 SUBROUTINE init_gphiv( ptab, i1, i2, j1, j2, before, nb,ndir) 453 use dom_oce 454 !!---------------------------------------------------------------------- 455 !! *** ROUTINE interpsshn *** 456 !!---------------------------------------------------------------------- 457 INTEGER , INTENT(in ) :: i1, i2, j1, j2 458 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 459 LOGICAL , INTENT(in ) :: before 460 INTEGER , INTENT(in ) :: nb , ndir 461 LOGICAL :: western_side, eastern_side,northern_side,southern_side 462 ! 463 !!---------------------------------------------------------------------- 464 ! 465 western_side = (nb == 1).AND.(ndir == 1) 466 eastern_side = (nb == 1).AND.(ndir == 2) 467 southern_side = (nb == 2).AND.(ndir == 1) 468 northern_side = (nb == 2).AND.(ndir == 2) 469 IF( before) THEN 470 ptab(i1:i2,j1:j2) = gphiv(i1:i2,j1:j2) 471 ELSE 472 gphiv(i1:i2,j1:j2)=ptab 473 ENDIF 474 ! 475 END SUBROUTINE init_gphiv 476 477 478 SUBROUTINE init_gphif( ptab, i1, i2, j1, j2, before, nb,ndir) 479 use dom_oce 480 !!---------------------------------------------------------------------- 481 !! *** ROUTINE interpsshn *** 482 !!---------------------------------------------------------------------- 483 INTEGER , INTENT(in ) :: i1, i2, j1, j2 484 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 485 LOGICAL , INTENT(in ) :: before 486 INTEGER , INTENT(in ) :: nb , ndir 487 LOGICAL :: western_side, eastern_side,northern_side,southern_side 488 ! 489 !!---------------------------------------------------------------------- 490 ! 491 western_side = (nb == 1).AND.(ndir == 1) 492 eastern_side = (nb == 1).AND.(ndir == 2) 493 southern_side = (nb == 2).AND.(ndir == 1) 494 northern_side = (nb == 2).AND.(ndir == 2) 495 IF( before) THEN 496 ptab(i1:i2,j1:j2) = gphif(i1:i2,j1:j2) 497 ELSE 498 gphif(i1:i2,j1:j2)=ptab 499 ENDIF 500 ! 501 END SUBROUTINE init_gphif 502 503 504 SUBROUTINE init_e1t( ptab, i1, i2, j1, j2, before, nb,ndir) 505 use dom_oce 506 !!---------------------------------------------------------------------- 507 !! *** ROUTINE interpsshn *** 508 !!---------------------------------------------------------------------- 509 INTEGER , INTENT(in ) :: i1, i2, j1, j2 510 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 511 LOGICAL , INTENT(in ) :: before 512 INTEGER , INTENT(in ) :: nb , ndir 513 LOGICAL :: western_side, eastern_side,northern_side,southern_side 514 ! 515 !!---------------------------------------------------------------------- 516 ! 517 western_side = (nb == 1).AND.(ndir == 1) 518 eastern_side = (nb == 1).AND.(ndir == 2) 519 southern_side = (nb == 2).AND.(ndir == 1) 520 northern_side = (nb == 2).AND.(ndir == 2) 521 IF( before) THEN 522 ptab(i1:i2,j1:j2) = e1t(i1:i2,j1:j2) 523 ELSE 524 e1t(i1:i2,j1:j2)=ptab/Agrif_rhoy() 525 ENDIF 526 ! 527 END SUBROUTINE init_e1t 528 529 SUBROUTINE init_e1u( ptab, i1, i2, j1, j2, before, nb,ndir) 530 use dom_oce 531 !!---------------------------------------------------------------------- 532 !! *** ROUTINE interpsshn *** 533 !!---------------------------------------------------------------------- 534 INTEGER , INTENT(in ) :: i1, i2, j1, j2 535 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 536 LOGICAL , INTENT(in ) :: before 537 INTEGER , INTENT(in ) :: nb , ndir 538 LOGICAL :: western_side, eastern_side,northern_side,southern_side 539 ! 540 !!---------------------------------------------------------------------- 541 ! 542 western_side = (nb == 1).AND.(ndir == 1) 543 eastern_side = (nb == 1).AND.(ndir == 2) 544 southern_side = (nb == 2).AND.(ndir == 1) 545 northern_side = (nb == 2).AND.(ndir == 2) 546 IF( before) THEN 547 ptab(i1:i2,j1:j2) = e1u(i1:i2,j1:j2) 548 ELSE 549 e1u(i1:i2,j1:j2)=ptab/Agrif_rhoy() 550 ENDIF 551 ! 552 END SUBROUTINE init_e1u 553 554 SUBROUTINE init_e1v( ptab, i1, i2, j1, j2, before, nb,ndir) 555 use dom_oce 556 !!---------------------------------------------------------------------- 557 !! *** ROUTINE interpsshn *** 558 !!---------------------------------------------------------------------- 559 INTEGER , INTENT(in ) :: i1, i2, j1, j2 560 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 561 LOGICAL , INTENT(in ) :: before 562 INTEGER , INTENT(in ) :: nb , ndir 563 LOGICAL :: western_side, eastern_side,northern_side,southern_side 564 ! 565 !!---------------------------------------------------------------------- 566 ! 567 western_side = (nb == 1).AND.(ndir == 1) 568 eastern_side = (nb == 1).AND.(ndir == 2) 569 southern_side = (nb == 2).AND.(ndir == 1) 570 northern_side = (nb == 2).AND.(ndir == 2) 571 IF( before) THEN 572 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) 573 ELSE 574 e1v(i1:i2,j1:j2)=ptab/Agrif_rhoy() 575 ENDIF 576 ! 577 END SUBROUTINE init_e1v 578 579 SUBROUTINE init_e1f( ptab, i1, i2, j1, j2, before, nb,ndir) 580 use dom_oce 581 !!---------------------------------------------------------------------- 582 !! *** ROUTINE interpsshn *** 583 !!---------------------------------------------------------------------- 584 INTEGER , INTENT(in ) :: i1, i2, j1, j2 585 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 586 LOGICAL , INTENT(in ) :: before 587 INTEGER , INTENT(in ) :: nb , ndir 588 LOGICAL :: western_side, eastern_side,northern_side,southern_side 589 ! 590 !!---------------------------------------------------------------------- 591 ! 592 western_side = (nb == 1).AND.(ndir == 1) 593 eastern_side = (nb == 1).AND.(ndir == 2) 594 southern_side = (nb == 2).AND.(ndir == 1) 595 northern_side = (nb == 2).AND.(ndir == 2) 596 IF( before) THEN 597 ptab(i1:i2,j1:j2) = e1f(i1:i2,j1:j2) 598 ELSE 599 e1f(i1:i2,j1:j2)=ptab/Agrif_rhoy() 600 ENDIF 601 ! 602 END SUBROUTINE init_e1f 603 604 SUBROUTINE init_e2t( ptab, i1, i2, j1, j2, before, nb,ndir) 605 use dom_oce 606 !!---------------------------------------------------------------------- 607 !! *** ROUTINE interpsshn *** 608 !!---------------------------------------------------------------------- 609 INTEGER , INTENT(in ) :: i1, i2, j1, j2 610 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 611 LOGICAL , INTENT(in ) :: before 612 INTEGER , INTENT(in ) :: nb , ndir 613 LOGICAL :: western_side, eastern_side,northern_side,southern_side 614 ! 615 !!---------------------------------------------------------------------- 616 ! 617 western_side = (nb == 1).AND.(ndir == 1) 618 eastern_side = (nb == 1).AND.(ndir == 2) 619 southern_side = (nb == 2).AND.(ndir == 1) 620 northern_side = (nb == 2).AND.(ndir == 2) 621 IF( before) THEN 622 ptab(i1:i2,j1:j2) = e2t(i1:i2,j1:j2) 623 ELSE 624 e2t(i1:i2,j1:j2)=ptab/Agrif_rhoy() 625 ENDIF 626 ! 627 END SUBROUTINE init_e2t 628 629 SUBROUTINE init_e2u( ptab, i1, i2, j1, j2, before, nb,ndir) 630 use dom_oce 631 !!---------------------------------------------------------------------- 632 !! *** ROUTINE interpsshn *** 633 !!---------------------------------------------------------------------- 634 INTEGER , INTENT(in ) :: i1, i2, j1, j2 635 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 636 LOGICAL , INTENT(in ) :: before 637 INTEGER , INTENT(in ) :: nb , ndir 638 LOGICAL :: western_side, eastern_side,northern_side,southern_side 639 ! 640 !!---------------------------------------------------------------------- 641 ! 642 western_side = (nb == 1).AND.(ndir == 1) 643 eastern_side = (nb == 1).AND.(ndir == 2) 644 southern_side = (nb == 2).AND.(ndir == 1) 645 northern_side = (nb == 2).AND.(ndir == 2) 646 IF( before) THEN 647 ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) 648 ELSE 649 e2u(i1:i2,j1:j2)=ptab/Agrif_rhoy() 650 ENDIF 651 ! 652 END SUBROUTINE init_e2u 653 654 SUBROUTINE init_e2v( ptab, i1, i2, j1, j2, before, nb,ndir) 655 use dom_oce 656 !!---------------------------------------------------------------------- 657 !! *** ROUTINE interpsshn *** 658 !!---------------------------------------------------------------------- 659 INTEGER , INTENT(in ) :: i1, i2, j1, j2 660 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 661 LOGICAL , INTENT(in ) :: before 662 INTEGER , INTENT(in ) :: nb , ndir 663 LOGICAL :: western_side, eastern_side,northern_side,southern_side 664 ! 665 !!---------------------------------------------------------------------- 666 ! 667 western_side = (nb == 1).AND.(ndir == 1) 668 eastern_side = (nb == 1).AND.(ndir == 2) 669 southern_side = (nb == 2).AND.(ndir == 1) 670 northern_side = (nb == 2).AND.(ndir == 2) 671 IF( before) THEN 672 ptab(i1:i2,j1:j2) = e2v(i1:i2,j1:j2) 673 ELSE 674 e2v(i1:i2,j1:j2)=ptab/Agrif_rhoy() 675 ENDIF 676 ! 677 END SUBROUTINE init_e2v 678 679 SUBROUTINE init_e2f( ptab, i1, i2, j1, j2, before, nb,ndir) 680 use dom_oce 681 !!---------------------------------------------------------------------- 682 !! *** ROUTINE interpsshn *** 683 !!---------------------------------------------------------------------- 684 INTEGER , INTENT(in ) :: i1, i2, j1, j2 685 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 686 LOGICAL , INTENT(in ) :: before 687 INTEGER , INTENT(in ) :: nb , ndir 688 LOGICAL :: western_side, eastern_side,northern_side,southern_side 689 ! 690 !!---------------------------------------------------------------------- 691 ! 692 western_side = (nb == 1).AND.(ndir == 1) 693 eastern_side = (nb == 1).AND.(ndir == 2) 694 southern_side = (nb == 2).AND.(ndir == 1) 695 northern_side = (nb == 2).AND.(ndir == 2) 1044 !!---------------------------------------------------------------------- 1045 INTEGER , INTENT(in ) :: i1, i2, j1, j2 1046 REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1047 LOGICAL , INTENT(in ) :: before 1048 INTEGER , INTENT(in ) :: nb , ndir 1049 ! 1050 !!---------------------------------------------------------------------- 1051 ! 696 1052 IF( before) THEN 697 1053 ptab(i1:i2,j1:j2) = e2f(i1:i2,j1:j2) … … 703 1059 704 1060 705 SUBROUTINE agrif_nemo_init706 USE agrif_parameters707 USE in_out_manager 708 USE lib_mpp 709 710 711 !!712 IMPLICIT NONE 713 714 INTEGER :: ios 715 716 NAMELIST/namagrif/ nn_cln_update,ln_spc_dyn,rn_sponge_tra,rn_sponge_dyn,ln_chk_bathy,npt_connect, npt_copy1061 SUBROUTINE agrif_nemo_init 1062 USE agrif_parameters 1063 USE dom_oce 1064 USE in_out_manager 1065 USE lib_mpp 1066 !! 1067 IMPLICIT NONE 1068 1069 INTEGER :: ios 1070 1071 NAMELIST/namagrif/ nn_cln_update,ln_spc_dyn,rn_sponge_tra,rn_sponge_dyn,ln_chk_bathy,npt_connect, & 1072 & npt_copy, ln_bry_south 717 1073 718 1074 REWIND( numnam_ref ) ! Namelist namagrif in reference namelist : nesting parameters … … 730 1086 WRITE(numout,*) '~~~~~~~' 731 1087 WRITE(numout,*) ' Namelist namagrif : set nesting parameters' 732 WRITE(numout,*) ' npt_copy = ', npt_copy 733 WRITE(numout,*) ' npt_connect = ', npt_connect 734 ENDIF 735 736 END SUBROUTINE agrif_nemo_init 737 738 739 SUBROUTINE Agrif_detect( kg, ksizex ) 1088 WRITE(numout,*) ' npt_copy = ', npt_copy 1089 WRITE(numout,*) ' npt_connect = ', npt_connect 1090 WRITE(numout,*) ' ln_bry_south = ', ln_bry_south 1091 ENDIF 1092 1093 ! Set the number of ghost cells according to periodicity 1094 1095 nbghostcells_x = nbghostcells 1096 nbghostcells_y_s = nbghostcells 1097 nbghostcells_y_n = nbghostcells 1098 1099 IF (.not.agrif_root()) THEN 1100 IF (jperio == 1) THEN 1101 nbghostcells_x = 0 1102 ENDIF 1103 IF (.NOT.ln_bry_south) THEN 1104 nbghostcells_y_s = 0 1105 ENDIF 1106 ENDIF 1107 1108 END SUBROUTINE agrif_nemo_init 1109 1110 SUBROUTINE Agrif_detect( kg, ksizex ) 740 1111 !!---------------------------------------------------------------------- 741 1112 !! *** ROUTINE Agrif_detect *** 742 1113 !!---------------------------------------------------------------------- 743 INTEGER, DIMENSION(2) :: ksizex 744 INTEGER, DIMENSION(ksizex(1),ksizex(2)) :: kg 745 !!---------------------------------------------------------------------- 746 ! 747 RETURN 748 ! 749 END SUBROUTINE Agrif_detect 750 SUBROUTINE agrif_before_regridding 751 END SUBROUTINE agrif_before_regridding 752 753 SUBROUTINE Agrif_InvLoc( indloc, nprocloc, i, indglob ) 754 !!---------------------------------------------------------------------- 755 !! *** ROUTINE Agrif_InvLoc *** 756 !!---------------------------------------------------------------------- 757 USE dom_oce 758 !! 759 IMPLICIT NONE 760 ! 761 INTEGER :: indglob, indloc, nprocloc, i 762 !!---------------------------------------------------------------------- 763 ! 764 SELECT CASE( i ) 765 CASE(1) ; indglob = indloc + nimppt(nprocloc+1) - 1 766 CASE(2) ; indglob = indloc + njmppt(nprocloc+1) - 1 767 CASE DEFAULT 768 indglob = indloc 769 END SELECT 770 ! 771 END SUBROUTINE Agrif_InvLoc 772 773 SUBROUTINE Agrif_get_proc_info( imin, imax, jmin, jmax ) 1114 INTEGER, DIMENSION(2) :: ksizex 1115 INTEGER, DIMENSION(ksizex(1),ksizex(2)) :: kg 1116 !!---------------------------------------------------------------------- 1117 ! 1118 RETURN 1119 ! 1120 END SUBROUTINE Agrif_detect 1121 1122 SUBROUTINE agrif_before_regridding 1123 END SUBROUTINE agrif_before_regridding 1124 1125 # if defined key_mpp_mpi 1126 SUBROUTINE Agrif_InvLoc( indloc, nprocloc, i, indglob ) 1127 !!---------------------------------------------------------------------- 1128 !! *** ROUTINE Agrif_InvLoc *** 1129 !!---------------------------------------------------------------------- 1130 USE dom_oce 1131 !! 1132 IMPLICIT NONE 1133 ! 1134 INTEGER :: indglob, indloc, nprocloc, i 1135 !!---------------------------------------------------------------------- 1136 ! 1137 SELECT CASE( i ) 1138 CASE(1) ; indglob = indloc + nimppt(nprocloc+1) - 1 1139 CASE(2) ; indglob = indloc + njmppt(nprocloc+1) - 1 1140 CASE DEFAULT 1141 indglob = indloc 1142 END SELECT 1143 ! 1144 END SUBROUTINE Agrif_InvLoc 1145 1146 SUBROUTINE Agrif_get_proc_info( imin, imax, jmin, jmax ) 774 1147 !!---------------------------------------------------------------------- 775 1148 !! *** ROUTINE Agrif_get_proc_info *** 776 1149 !!---------------------------------------------------------------------- 777 USE par_oce778 USE dom_oce779 !!780 IMPLICIT NONE781 !782 INTEGER, INTENT(out) :: imin, imax783 INTEGER, INTENT(out) :: jmin, jmax784 !!----------------------------------------------------------------------785 !786 imin = nimppt(Agrif_Procrank+1) ! ?????787 jmin = njmppt(Agrif_Procrank+1) ! ?????788 imax = imin + jpi - 1789 jmax = jmin + jpj - 1790 !791 END SUBROUTINE Agrif_get_proc_info792 793 SUBROUTINE Agrif_estimate_parallel_cost(imin, imax,jmin, jmax, nbprocs, grid_cost)1150 USE par_oce 1151 USE dom_oce 1152 !! 1153 IMPLICIT NONE 1154 ! 1155 INTEGER, INTENT(out) :: imin, imax 1156 INTEGER, INTENT(out) :: jmin, jmax 1157 !!---------------------------------------------------------------------- 1158 ! 1159 imin = nimppt(Agrif_Procrank+1) ! ????? 1160 jmin = njmppt(Agrif_Procrank+1) ! ????? 1161 imax = imin + jpi - 1 1162 jmax = jmin + jpj - 1 1163 ! 1164 END SUBROUTINE Agrif_get_proc_info 1165 1166 SUBROUTINE Agrif_estimate_parallel_cost(imin, imax,jmin, jmax, nbprocs, grid_cost) 794 1167 !!---------------------------------------------------------------------- 795 1168 !! *** ROUTINE Agrif_estimate_parallel_cost *** 796 1169 !!---------------------------------------------------------------------- 797 USE par_oce798 !!799 IMPLICIT NONE800 !801 INTEGER, INTENT(in) :: imin, imax802 INTEGER, INTENT(in) :: jmin, jmax803 INTEGER, INTENT(in) :: nbprocs804 REAL(wp), INTENT(out) :: grid_cost805 !!----------------------------------------------------------------------806 !807 grid_cost = REAL(imax-imin+1,wp)*REAL(jmax-jmin+1,wp) / REAL(nbprocs,wp)808 !809 END SUBROUTINE Agrif_estimate_parallel_cost810 1170 USE par_oce 1171 !! 1172 IMPLICIT NONE 1173 ! 1174 INTEGER, INTENT(in) :: imin, imax 1175 INTEGER, INTENT(in) :: jmin, jmax 1176 INTEGER, INTENT(in) :: nbprocs 1177 REAL(wp), INTENT(out) :: grid_cost 1178 !!---------------------------------------------------------------------- 1179 ! 1180 grid_cost = REAL(imax-imin+1,wp)*REAL(jmax-jmin+1,wp) / REAL(nbprocs,wp) 1181 ! 1182 END SUBROUTINE Agrif_estimate_parallel_cost 1183 # endif 811 1184 #endif -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/bilinear_interp.f90
r12414 r13024 1 ! 1 2 MODULE bilinear_interp 2 3 END MODULE 3 ! 4 USE agrif_modutil 5 ! 6 !************************************************************************ 7 ! * 8 ! MODULE BILINEAR INTERP * 9 ! * 10 ! bilinear interpolation routines from SCRIP package * 11 ! * 12 !http://climate.lanl.gov/Software/SCRIP/ * 13 ! * 14 !Bilinear remapping * 15 ! * 16 !************************************************************************ 17 ! 18 !----------------------------------------------------------------------- 19 IMPLICIT NONE 20 21 !----------------------------------------------------------------------- 22 ! variables that describe each grid 23 !----------------------------------------------------------------------- 24 ! 25 INTEGER :: grid1_size,grid2_size,grid1_rank, grid2_rank 26 ! 27 INTEGER, DIMENSION(:), ALLOCATABLE :: grid1_dims, grid2_dims 28 ! 29 30 !----------------------------------------------------------------------- 31 ! grid coordinates and masks 32 !----------------------------------------------------------------------- 33 ! 34 LOGICAL, DIMENSION(:), ALLOCATABLE :: grid1_mask,grid2_mask 35 ! each grid center in radians 36 REAL*8,DIMENSION(:), ALLOCATABLE :: & 37 grid1_center_lat, & 38 grid1_center_lon, & 39 grid2_center_lat, & 40 grid2_center_lon, & 41 grid1_frac, & ! fractional area of grid cells 42 grid2_frac ! participating in remapping 43 ! 44 ! lat/lon bounding box for use in restricting grid searches 45 ! 46 REAL*8,DIMENSION(:,:), ALLOCATABLE :: grid1_bound_box,grid2_bound_box 47 ! 48 !----------------------------------------------------------------------- 49 ! bins for restricting searches 50 !----------------------------------------------------------------------- 51 ! 52 ! num of bins for restricted srch 53 INTEGER, PARAMETER :: num_srch_bins = 90 54 ! 55 ! min,max adds for grid cells in this lat bin 56 ! 57 INTEGER,DIMENSION(:,:), ALLOCATABLE :: bin_addr1,bin_addr2 58 ! 59 ! min,max longitude for each search bin 60 ! 61 REAL*8, DIMENSION(:,:), ALLOCATABLE :: bin_lats,bin_lons 62 63 REAL*8, PARAMETER :: zero = 0.0, & 64 one = 1.0, & 65 two = 2.0, & 66 three = 3.0, & 67 four = 4.0, & 68 five = 5.0, & 69 half = 0.5, & 70 quart = 0.25, & 71 bignum = 1.e+20, & 72 tiny = 1.e-14, & 73 pi = 3.14159265359, & 74 pi2 = two*pi, & 75 pih = half*pi 76 77 REAL*8, PARAMETER :: deg2rad = pi/180. 78 ! 79 ! max iteration count for i,j iteration 80 ! 81 INTEGER , PARAMETER :: max_iter = 100 82 ! 83 ! convergence criterion 84 ! 85 REAL*8, PARAMETER :: converge = 1.e-10 86 87 88 ! 89 INTEGER, PARAMETER :: norm_opt_none = 1 & 90 ,norm_opt_dstarea = 2 & 91 ,norm_opt_frcarea = 3 92 ! 93 INTEGER, PARAMETER :: map_type_conserv = 1 & 94 ,map_type_bilinear = 2 & 95 ,map_type_bicubic = 3 & 96 ,map_type_distwgt = 4 97 ! 98 INTEGER :: max_links_map1 & ! current size of link arrays 99 ,num_links_map1 & ! actual number of links for remapping 100 ,max_links_map2 & ! current size of link arrays 101 ,num_links_map2 & ! actual number of links for remapping 102 ,num_maps & ! num of remappings for this grid pair 103 ,num_wts & ! num of weights used in remapping 104 ,map_type & ! identifier for remapping method 105 ,norm_opt & ! option for normalization (conserv only) 106 ,resize_increment ! default amount to increase array size 107 108 INTEGER , DIMENSION(:), ALLOCATABLE :: & 109 grid1_add_map1, & ! grid1 address for each link in mapping 1 110 grid2_add_map1, & ! grid2 address for each link in mapping 1 111 grid1_add_map2, & ! grid1 address for each link in mapping 2 112 grid2_add_map2 ! grid2 address for each link in mapping 2 113 114 REAL*8, DIMENSION(:,:), ALLOCATABLE :: & 115 wts_map1, & ! map weights for each link (num_wts,max_links) 116 wts_map2 ! map weights for each link (num_wts,max_links) 117 ! 118 CONTAINS 119 ! 120 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 121 !************************************************************************ 122 ! SUBROUTINE GRID_INIT 123 !************************************************************************ 124 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 125 ! 126 SUBROUTINE get_remap_matrix(grid1_lat,grid2_lat,grid1_lon,grid2_lon,mask, & 127 remap_matrix,source_add,destination_add) 128 ! 129 !----------------------------------------------------------------------- 130 !this routine makes any necessary changes (e.g. for 0,2pi longitude range) 131 !----------------------------------------------------------------------- 132 ! 133 REAL*8,DIMENSION(:,:) :: grid1_lat,grid2_lat,grid1_lon,grid2_lon 134 LOGICAL,DIMENSION(:,:) :: mask 135 ! 136 INTEGER,DIMENSION(:),ALLOCATABLE :: source_add,destination_add 137 REAL*8,DIMENSION(:,:), ALLOCATABLE :: remap_matrix 138 ! 139 !----------------------------------------------------------------------- 140 ! local variables 141 !----------------------------------------------------------------------- 142 ! 143 INTEGER :: n,nele,i,j,ip1,jp1,n_add,e_add,ne_add,nx,ny 144 INTEGER :: xpos,ypos 145 ! 146 ! integer mask 147 ! 148 INTEGER, DIMENSION(:), ALLOCATABLE :: imask 149 ! 150 ! lat/lon intervals for search bins 151 ! 152 REAL*8 :: dlat,dlon 153 ! 154 ! temps for computing bounding boxes 155 ! 156 REAL*8, DIMENSION(4) :: tmp_lats, tmp_lons 157 ! 158 ! write(*,*)'proceed to Bilinear interpolation ...' 159 ! 160 ! IF(ALLOCATED(wts_map1)) DEALLOCATE(wts_map1) 161 ! IF(ALLOCATED(grid1_add_map1)) DEALLOCATE(grid1_add_map1) 162 ! IF(ALLOCATED(grid2_add_map1)) DEALLOCATE(grid2_add_map1) 163 164 165 ! 166 ALLOCATE(grid1_dims(2),grid2_dims(2)) 167 ! 168 grid1_dims(1) = SIZE(grid1_lat,2) 169 grid1_dims(2) = SIZE(grid1_lat,1) 170 grid2_dims(1) = SIZE(grid2_lat,2) 171 grid2_dims(2) = SIZE(grid2_lat,1) 172 grid1_size = SIZE(grid1_lat,2) * SIZE(grid1_lat,1) 173 grid2_size = SIZE(grid2_lat,2) * SIZE(grid2_lat,1) 174 ! 175 !----------------------------------------------------------------------- 176 ! allocate grid coordinates/masks and read data 177 !----------------------------------------------------------------------- 178 ! 179 ALLOCATE( grid2_mask(grid2_size), & 180 grid1_bound_box (4,grid1_size), & 181 grid2_bound_box (4,grid2_size), & 182 grid1_frac (grid1_size), & 183 grid2_frac (grid2_size)) 184 ALLOCATE(imask(grid1_size)) 185 ! 186 ! 187 grid1_frac = zero 188 grid2_frac = zero 189 190 ! 191 ! 2D array -> 1D array 192 ! 193 ALLOCATE(grid1_center_lat(SIZE(grid1_lat,1)*SIZE(grid1_lat,2))) 194 CALL tab2Dto1D(grid1_lat,grid1_center_lat) 195 196 ALLOCATE(grid1_center_lon(SIZE(grid1_lon,1)*SIZE(grid1_lon,2))) 197 CALL tab2Dto1D(grid1_lon,grid1_center_lon) 198 199 ALLOCATE(grid2_center_lat(SIZE(grid2_lat,1)*SIZE(grid2_lat,2))) 200 CALL tab2Dto1D(grid2_lat,grid2_center_lat) 201 202 ALLOCATE(grid2_center_lon(SIZE(grid2_lon,1)*SIZE(grid2_lon,2))) 203 CALL tab2Dto1D(grid2_lon,grid2_center_lon) 204 205 ALLOCATE(grid1_mask(SIZE(grid1_lat,1)*SIZE(grid1_lat,2))) 206 CALL logtab2Dto1D(mask,grid1_mask) 207 ! 208 ! Write(*,*) ,'grid1_mask = ',grid1_mask 209 ! 210 ! degrees to radian 211 ! 212 grid1_center_lat = grid1_center_lat*deg2rad 213 grid1_center_lon = grid1_center_lon*deg2rad 214 grid2_center_lat = grid2_center_lat*deg2rad 215 grid2_center_lon = grid2_center_lon*deg2rad 216 217 !----------------------------------------------------------------------- 218 ! convert longitudes to 0,2pi interval 219 !----------------------------------------------------------------------- 220 221 WHERE (grid1_center_lon .GT. pi2) grid1_center_lon = & 222 grid1_center_lon - pi2 223 WHERE (grid1_center_lon .LT. zero) grid1_center_lon = & 224 grid1_center_lon + pi2 225 WHERE (grid2_center_lon .GT. pi2) grid2_center_lon = & 226 grid2_center_lon - pi2 227 WHERE (grid2_center_lon .LT. zero) grid2_center_lon = & 228 grid2_center_lon + pi2 229 230 !----------------------------------------------------------------------- 231 ! 232 ! make sure input latitude range is within the machine values 233 ! for +/- pi/2 234 ! 235 !----------------------------------------------------------------------- 236 237 WHERE (grid1_center_lat > pih) grid1_center_lat = pih 238 WHERE (grid1_center_lat < -pih) grid1_center_lat = -pih 239 WHERE (grid2_center_lat > pih) grid2_center_lat = pih 240 WHERE (grid2_center_lat < -pih) grid2_center_lat = -pih 241 242 !----------------------------------------------------------------------- 243 ! 244 ! compute bounding boxes for restricting future grid searches 245 ! 246 !----------------------------------------------------------------------- 247 ! 248 nx = grid1_dims(1) 249 ny = grid1_dims(2) 250 251 DO n=1,grid1_size 252 253 !*** find N,S and NE points to this grid point 254 255 j = (n - 1)/nx +1 256 i = n - (j-1)*nx 257 258 IF (i < nx) THEN 259 ip1 = i + 1 260 ELSE 261 !*** assume cyclic 262 ip1 = 1 263 !*** but if it is not, correct 264 e_add = (j - 1)*nx + ip1 265 IF (ABS(grid1_center_lat(e_add) - & 266 grid1_center_lat(n )) > pih) THEN 267 ip1 = i 268 ENDIF 269 ip1=nx 270 ENDIF 271 272 IF (j < ny) THEN 273 jp1 = j+1 274 ELSE 275 !*** assume cyclic 276 jp1 = 1 277 !*** but if it is not, correct 278 n_add = (jp1 - 1)*nx + i 279 IF (ABS(grid1_center_lat(n_add) - & 280 grid1_center_lat(n )) > pih) THEN 281 jp1 = j 282 ENDIF 283 jp1=ny 284 ENDIF 285 286 n_add = (jp1 - 1)*nx + i 287 e_add = (j - 1)*nx + ip1 288 ne_add = (jp1 - 1)*nx + ip1 289 290 !*** find N,S and NE lat/lon coords and check bounding box 291 292 tmp_lats(1) = grid1_center_lat(n) 293 tmp_lats(2) = grid1_center_lat(e_add) 294 tmp_lats(3) = grid1_center_lat(ne_add) 295 tmp_lats(4) = grid1_center_lat(n_add) 296 297 tmp_lons(1) = grid1_center_lon(n) 298 tmp_lons(2) = grid1_center_lon(e_add) 299 tmp_lons(3) = grid1_center_lon(ne_add) 300 tmp_lons(4) = grid1_center_lon(n_add) 301 302 grid1_bound_box(1,n) = MINVAL(tmp_lats) 303 grid1_bound_box(2,n) = MAXVAL(tmp_lats) 304 305 grid1_bound_box(3,n) = MINVAL(tmp_lons) 306 grid1_bound_box(4,n) = MAXVAL(tmp_lons) 307 END DO 308 309 nx = grid2_dims(1) 310 ny = grid2_dims(2) 311 312 DO n=1,grid2_size 313 314 !*** find N,S and NE points to this grid point 315 316 j = (n - 1)/nx +1 317 i = n - (j-1)*nx 318 319 IF (i < nx) THEN 320 ip1 = i + 1 321 ELSE 322 !*** assume cyclic 323 ip1 = 1 324 !*** but if it is not, correct 325 e_add = (j - 1)*nx + ip1 326 IF (ABS(grid2_center_lat(e_add) - & 327 grid2_center_lat(n )) > pih) THEN 328 ip1 = i 329 ENDIF 330 ENDIF 331 332 IF (j < ny) THEN 333 jp1 = j+1 334 ELSE 335 !*** assume cyclic 336 jp1 = 1 337 !*** but if it is not, correct 338 n_add = (jp1 - 1)*nx + i 339 IF (ABS(grid2_center_lat(n_add) - & 340 grid2_center_lat(n )) > pih) THEN 341 jp1 = j 342 ENDIF 343 ENDIF 344 345 n_add = (jp1 - 1)*nx + i 346 e_add = (j - 1)*nx + ip1 347 ne_add = (jp1 - 1)*nx + ip1 348 349 !*** find N,S and NE lat/lon coords and check bounding box 350 351 tmp_lats(1) = grid2_center_lat(n) 352 tmp_lats(2) = grid2_center_lat(e_add) 353 tmp_lats(3) = grid2_center_lat(ne_add) 354 tmp_lats(4) = grid2_center_lat(n_add) 355 356 tmp_lons(1) = grid2_center_lon(n) 357 tmp_lons(2) = grid2_center_lon(e_add) 358 tmp_lons(3) = grid2_center_lon(ne_add) 359 tmp_lons(4) = grid2_center_lon(n_add) 360 361 grid2_bound_box(1,n) = MINVAL(tmp_lats) 362 grid2_bound_box(2,n) = MAXVAL(tmp_lats) 363 grid2_bound_box(3,n) = MINVAL(tmp_lons) 364 grid2_bound_box(4,n) = MAXVAL(tmp_lons) 365 END DO 366 ! 367 ! 368 ! 369 WHERE (ABS(grid1_bound_box(4,:) - grid1_bound_box(3,:)) > pi) 370 grid1_bound_box(3,:) = zero 371 grid1_bound_box(4,:) = pi2 372 END WHERE 373 374 WHERE (ABS(grid2_bound_box(4,:) - grid2_bound_box(3,:)) > pi) 375 grid2_bound_box(3,:) = zero 376 grid2_bound_box(4,:) = pi2 377 END WHERE 378 379 !*** 380 !*** try to check for cells that overlap poles 381 !*** 382 383 WHERE (grid1_center_lat > grid1_bound_box(2,:)) & 384 grid1_bound_box(2,:) = pih 385 386 WHERE (grid1_center_lat < grid1_bound_box(1,:)) & 387 grid1_bound_box(1,:) = -pih 388 389 WHERE (grid2_center_lat > grid2_bound_box(2,:)) & 390 grid2_bound_box(2,:) = pih 391 392 WHERE (grid2_center_lat < grid2_bound_box(1,:)) & 393 grid2_bound_box(1,:) = -pih 394 395 !----------------------------------------------------------------------- 396 ! set up and assign address ranges to search bins in order to 397 ! further restrict later searches 398 !----------------------------------------------------------------------- 399 400 ALLOCATE(bin_addr1(2,num_srch_bins)) 401 ALLOCATE(bin_addr2(2,num_srch_bins)) 402 ALLOCATE(bin_lats (2,num_srch_bins)) 403 ALLOCATE(bin_lons (2,num_srch_bins)) 404 405 dlat = pi/num_srch_bins 406 407 DO n=1,num_srch_bins 408 bin_lats(1,n) = (n-1)*dlat - pih 409 bin_lats(2,n) = n*dlat - pih 410 bin_lons(1,n) = zero 411 bin_lons(2,n) = pi2 412 bin_addr1(1,n) = grid1_size + 1 413 bin_addr1(2,n) = 0 414 bin_addr2(1,n) = grid2_size + 1 415 bin_addr2(2,n) = 0 416 END DO 417 418 DO nele=1,grid1_size 419 DO n=1,num_srch_bins 420 IF (grid1_bound_box(1,nele) <= bin_lats(2,n) .AND. & 421 grid1_bound_box(2,nele) >= bin_lats(1,n)) THEN 422 bin_addr1(1,n) = MIN(nele,bin_addr1(1,n)) 423 bin_addr1(2,n) = MAX(nele,bin_addr1(2,n)) 424 ENDIF 425 END DO 426 END DO 427 428 DO nele=1,grid2_size 429 DO n=1,num_srch_bins 430 IF (grid2_bound_box(1,nele) <= bin_lats(2,n) .AND. & 431 grid2_bound_box(2,nele) >= bin_lats(1,n)) THEN 432 bin_addr2(1,n) = MIN(nele,bin_addr2(1,n)) 433 bin_addr2(2,n) = MAX(nele,bin_addr2(2,n)) 434 ENDIF 435 END DO 436 END DO 437 ! 438 CALL init_remap_vars 439 CALL remap_bilin 440 441 ALLOCATE(remap_matrix(SIZE(wts_map1,1),SIZE(wts_map1,2)), & 442 source_add(SIZE(grid1_add_map1)), & 443 destination_add(SIZE(grid2_add_map1))) 444 445 DO j = 1,SIZE(wts_map1,2) 446 DO i = 1,SIZE(wts_map1,1) 447 448 remap_matrix(i,j) = wts_map1(i,j) 449 450 END DO 451 END DO 452 453 454 source_add(:) = grid1_add_map1(:) 455 destination_add(:) = grid2_add_map1(:) 456 ! 457 WHERE(destination_add == 0) 458 destination_add = 1 459 END WHERE 460 461 WHERE(source_add == 0) 462 source_add = 1 463 END WHERE 464 ! 465 !DEALLOCATE(grid1_bound_box,grid2_bound_box,grid1_center_lat,grid1_center_lon) 466 !DEALLOCATE(grid2_center_lat,grid2_center_lon,grid2_add_map1,grid1_add_map1,wts_map1) 467 !DEALLOCATE(grid1_frac,grid2_frac,grid1_dims,grid2_dims,grid2_mask,imask) 468 !DEALLOCATE(bin_addr1,bin_addr2,bin_lats,bin_lons) 469 !DEALLOCATE(grid1_mask) 470 ! 471 !----------------------------------------------------------------------- 472 473 END SUBROUTINE get_remap_matrix 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 !*********************************************************************** 499 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 500 !************************************************************************ 501 ! SUBROUTINE REMAP_BILINEAR 502 !************************************************************************ 503 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 504 505 SUBROUTINE remap_bilin 506 507 !----------------------------------------------------------------------- 508 ! this routine computes the weights for a bilinear interpolation. 509 !----------------------------------------------------------------------- 510 511 !----------------------------------------------------------------------- 512 ! local variables 513 !----------------------------------------------------------------------- 514 515 INTEGER :: n,icount,dst_add,iter,nmap,nbmasked 516 ! 517 ! address for the four source points 518 ! 519 INTEGER, DIMENSION(4) :: src_add,involved_pts 520 INTEGER, DIMENSION(1) :: minlon 521 INTEGER, DIMENSION(1) :: minlat 522 REAL*8, DIMENSION(4) :: distx,disty 523 REAL*8 :: normalize 524 ! 525 ! latitudes longitudes of four bilinear corners 526 ! 527 REAL*8, DIMENSION(4) :: src_lats,src_lons 528 ! 529 ! bilinear weights for four corners 530 ! 531 REAL*8, DIMENSION(4) :: wgts 532 ! 533 REAL*8 :: & 534 plat, plon, & ! lat/lon coords of destination point 535 iguess, jguess, & ! current guess for bilinear coordinate 536 thguess, phguess, & ! current guess for lat/lon coordinate 537 deli, delj, & ! corrections to i,j 538 dth1, dth2, dth3, & ! some latitude differences 539 dph1, dph2, dph3, & ! some longitude differences 540 dthp, dphp, & ! difference between point and sw corner 541 mat1, mat2, mat3, mat4, & ! matrix elements 542 determinant, & ! matrix determinant 543 sum_wgts ! sum of weights for normalization 544 545 INTEGER lastsrc_add 546 REAL*8,DIMENSION(:),ALLOCATABLE :: cos_grid1_center_lat 547 REAL*8,DIMENSION(:),ALLOCATABLE :: cos_grid1_center_lon 548 REAL*8,DIMENSION(:),ALLOCATABLE :: sin_grid1_center_lat 549 REAL*8,DIMENSION(:),ALLOCATABLE :: sin_grid1_center_lon 550 INTEGER :: i 551 ! 552 grid2_mask = .TRUE. 553 ! 554 ! 555 nmap = 1 556 ! 557 !*** 558 !*** loop over destination grid 559 !*** 560 ! print*,'grid2_size =',grid2_size 561 ! 562 lastsrc_add=1 563 ! 564 Allocate(cos_grid1_center_lat(size(grid1_center_lat))) 565 Allocate(sin_grid1_center_lat(size(grid1_center_lat))) 566 Allocate(cos_grid1_center_lon(size(grid1_center_lon))) 567 Allocate(sin_grid1_center_lon(size(grid1_center_lon))) 568 569 do i=1,size(grid1_center_lat) 570 cos_grid1_center_lat(i)=cos(grid1_center_lat(i)) 571 sin_grid1_center_lat(i)=sin(grid1_center_lat(i)) 572 ENDDO 573 574 do i=1,size(grid1_center_lon) 575 cos_grid1_center_lon(i)=cos(grid1_center_lon(i)) 576 sin_grid1_center_lon(i)=sin(grid1_center_lon(i)) 577 ENDDO 578 579 grid_loop1: DO dst_add = 1, grid2_size 580 581 IF (.NOT. grid2_mask(dst_add)) CYCLE grid_loop1 582 ! 583 plat = grid2_center_lat(dst_add) 584 plon = grid2_center_lon(dst_add) 585 !*** 586 !*** find nearest square of grid points on source grid 587 !*** 588 CALL grid_search_bilin(src_add, src_lats, src_lons, & 589 plat, plon, grid1_dims, & 590 grid1_center_lat, cos_grid1_center_lat, sin_grid1_center_lat, & 591 grid1_center_lon, cos_grid1_center_lon, sin_grid1_center_lon, & 592 grid1_bound_box, bin_addr1, bin_addr2,lastsrc_add) 593 !*** 594 !*** check to see if points are land points 595 !*** 596 ! 597 IF (src_add(1) > 0) THEN 598 ! 599 DO n=1,4 600 ! if(.not. grid1_mask(src_add(n))) nbmasked = nbmasked + 1 601 IF(.NOT. grid1_mask(src_add(n))) src_add(1) = 0 602 END DO 603 ! 604 ENDIF 605 ! 606 ! 607 608 !*** 609 !*** if point found, find local i,j coordinates for weights 610 !*** 611 IF (src_add(1) > 0) THEN 612 grid2_frac(dst_add) = one 613 !*** 614 !*** iterate to find i,j for bilinear approximation 615 !*** 616 dth1 = src_lats(2) - src_lats(1) 617 dth2 = src_lats(4) - src_lats(1) 618 dth3 = src_lats(3) - src_lats(2) - dth2 619 620 dph1 = src_lons(2) - src_lons(1) 621 dph2 = src_lons(4) - src_lons(1) 622 dph3 = src_lons(3) - src_lons(2) 623 624 IF (dph1 > three*pih) dph1 = dph1 - pi2 625 IF (dph2 > three*pih) dph2 = dph2 - pi2 626 IF (dph3 > three*pih) dph3 = dph3 - pi2 627 IF (dph1 < -three*pih) dph1 = dph1 + pi2 628 IF (dph2 < -three*pih) dph2 = dph2 + pi2 629 IF (dph3 < -three*pih) dph3 = dph3 + pi2 630 631 dph3 = dph3 - dph2 632 633 iguess = half 634 jguess = half 635 636 iter_loop1: DO iter=1,max_iter 637 638 dthp = plat - src_lats(1) - dth1*iguess - & 639 dth2*jguess - dth3*iguess*jguess 640 dphp = plon - src_lons(1) 641 642 IF (dphp > three*pih) dphp = dphp - pi2 643 IF (dphp < -three*pih) dphp = dphp + pi2 644 645 dphp = dphp - dph1*iguess - dph2*jguess - & 646 dph3*iguess*jguess 647 648 mat1 = dth1 + dth3*jguess 649 mat2 = dth2 + dth3*iguess 650 mat3 = dph1 + dph3*jguess 651 mat4 = dph2 + dph3*iguess 652 653 determinant = mat1*mat4 - mat2*mat3 654 655 deli = (dthp*mat4 - mat2*dphp)/determinant 656 delj = (mat1*dphp - dthp*mat3)/determinant 657 658 IF (ABS(deli) < converge .AND. & 659 ABS(delj) < converge) EXIT iter_loop1 660 661 iguess = iguess + deli 662 jguess = jguess + delj 663 664 END DO iter_loop1 665 666 IF (iter <= max_iter) THEN 667 668 !*** 669 !*** successfully found i,j - compute weights 670 !*** 671 672 wgts(1) = (one-iguess)*(one-jguess) 673 wgts(2) = iguess*(one-jguess) 674 wgts(3) = iguess*jguess 675 wgts(4) = (one-iguess)*jguess 676 ! 677 ! 678 CALL store_link_bilin(dst_add, src_add, wgts, nmap) 679 680 ELSE 681 PRINT *,'Point coords: ',plat,plon 682 PRINT *,'Dest grid lats: ',src_lats 683 PRINT *,'Dest grid lons: ',src_lons 684 PRINT *,'Dest grid addresses: ',src_add 685 PRINT *,'Current i,j : ',iguess, jguess 686 STOP 'Iteration for i,j exceed max iteration count' 687 ENDIF 688 ! 689 !*** 690 !*** search for bilinear failed - use a distance-weighted 691 !*** average instead (this is typically near the pole) 692 !*** 693 ELSE IF (src_add(1) < 0) THEN 694 695 src_add = ABS(src_add) 696 icount = 0 697 DO n=1,4 698 ! 699 IF (grid1_mask(src_add(n))) THEN 700 icount = icount + 1 701 ELSE 702 src_lats(n) = zero 703 ENDIF 704 ! 705 END DO 706 707 IF (icount > 0) THEN 708 ! 709 !*** renormalize weights 710 ! 711 sum_wgts = SUM(src_lats) 712 wgts(1) = src_lats(1)/sum_wgts 713 wgts(2) = src_lats(2)/sum_wgts 714 wgts(3) = src_lats(3)/sum_wgts 715 wgts(4) = src_lats(4)/sum_wgts 716 ! 717 grid2_frac(dst_add) = one 718 CALL store_link_bilin(dst_add, src_add, wgts, nmap) 719 ENDIF 720 721 ! 722 ENDIF 723 END DO grid_loop1 724 ! 725 ! Call sort_add(grid2_add_map1, grid1_add_map1, wts_map1) 726 ! 727 ! 728 !----------------------------------------------------------------------- 729 730 deallocate(cos_grid1_center_lat) 731 deallocate(cos_grid1_center_lon) 732 deallocate(sin_grid1_center_lat) 733 deallocate(sin_grid1_center_lon) 734 END SUBROUTINE remap_bilin 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 !*********************************************************************** 754 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 755 !************************************************************************ 756 ! SUBROUTINE GRID_SEARCH_BILIN 757 !************************************************************************ 758 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 759 ! 760 SUBROUTINE grid_search_bilin(src_add, src_lats, src_lons, & 761 plat, plon, src_grid_dims, & 762 src_center_lat, cos_src_center_lat, sin_src_center_lat, & 763 src_center_lon, cos_src_center_lon, sin_src_center_lon,& 764 src_grid_bound_box, & 765 src_bin_add, dst_bin_add,lastsrc_add) 766 767 !----------------------------------------------------------------------- 768 ! 769 ! this routine finds the location of the search point plat, plon 770 ! in the source grid and returns the corners needed for a bilinear 771 ! interpolation. 772 ! 773 !----------------------------------------------------------------------- 774 775 !----------------------------------------------------------------------- 776 ! output variables 777 !----------------------------------------------------------------------- 778 ! 779 ! address of each corner point enclosing P 780 ! 781 INTEGER,DIMENSION(4) :: src_add 782 REAL*8,DIMENSION(4) :: src_lats,src_lons 783 ! 784 !----------------------------------------------------------------------- 785 ! input variables 786 !----------------------------------------------------------------------- 787 ! 788 ! latitude, longitude of the search point 789 ! 790 REAL*8, INTENT(in) :: plat,plon 791 ! 792 ! size of each src grid dimension 793 ! 794 INTEGER, DIMENSION(2), INTENT(in) :: src_grid_dims 795 ! 796 ! latitude, longitude of each src grid center 797 ! 798 REAL*8, DIMENSION(:), INTENT(in) :: src_center_lat,src_center_lon 799 REAL*8, DIMENSION(:), INTENT(in) :: cos_src_center_lat 800 REAL*8, DIMENSION(:), INTENT(in) :: sin_src_center_lat 801 REAL*8, DIMENSION(:), INTENT(in) :: cos_src_center_lon 802 REAL*8, DIMENSION(:), INTENT(in) :: sin_src_center_lon 803 ! 804 ! bound box for source grid 805 ! 806 REAL*8, DIMENSION(:,:), INTENT(in) :: src_grid_bound_box 807 ! 808 ! latitude bins for restricting searches 809 ! 810 INTEGER, DIMENSION(:,:), INTENT(in) ::src_bin_add,dst_bin_add 811 812 INTEGER,OPTIONAL :: lastsrc_add 813 INTEGER :: loopsrc,l1,l2 814 ! 815 !----------------------------------------------------------------------- 816 ! local variables 817 !----------------------------------------------------------------------- 818 ! 819 INTEGER :: n,next_n,srch_add,nx, ny,min_add, max_add, & 820 i, j, jp1, ip1, n_add, e_add, ne_add 821 822 823 REAL*8 :: vec1_lat, vec1_lon,vec2_lat, vec2_lon, cross_product, & 824 cross_product_last,coslat_dst, sinlat_dst, coslon_dst, & 825 sinlon_dst,dist_min, distance 826 827 !----------------------------------------------------------------------- 828 ! restrict search first using bins 829 !----------------------------------------------------------------------- 830 831 src_add = 0 832 833 min_add = SIZE(src_center_lat) 834 max_add = 1 835 DO n=1,num_srch_bins 836 IF (plat >= bin_lats(1,n) .AND. plat <= bin_lats(2,n) .AND. & 837 plon >= bin_lons(1,n) .AND. plon <= bin_lons(2,n)) THEN 838 min_add = MIN(min_add, src_bin_add(1,n)) 839 max_add = MAX(max_add, src_bin_add(2,n)) 840 ENDIF 841 END DO 842 843 !----------------------------------------------------------------------- 844 ! now perform a more detailed search 845 !----------------------------------------------------------------------- 846 847 nx = src_grid_dims(1) 848 ny = src_grid_dims(2) 849 850 loopsrc=0 851 DO WHILE (loopsrc <= max_add) 852 853 854 l1=MAX(min_add,lastsrc_add-loopsrc) 855 l2=MIN(max_add,lastsrc_add+loopsrc) 856 857 loopsrc = loopsrc+1 858 859 srch_loop: DO srch_add = l1,l2,MAX(l2-l1,1) 860 861 !*** first check bounding box 862 863 IF (plat <= src_grid_bound_box(2,srch_add) .AND. & 864 plat >= src_grid_bound_box(1,srch_add) .AND. & 865 plon <= src_grid_bound_box(4,srch_add) .AND. & 866 plon >= src_grid_bound_box(3,srch_add)) THEN 867 !*** 868 !*** we are within bounding box so get really serious 869 !*** 870 !*** determine neighbor addresses 871 ! 872 j = (srch_add - 1)/nx +1 873 i = srch_add - (j-1)*nx 874 ! 875 IF (i < nx) THEN 876 ip1 = i + 1 877 ELSE 878 ip1 = 1 879 ENDIF 880 ! 881 IF (j < ny) THEN 882 jp1 = j+1 883 ELSE 884 jp1 = 1 885 ENDIF 886 ! 887 n_add = (jp1 - 1)*nx + i 888 e_add = (j - 1)*nx + ip1 889 ne_add = (jp1 - 1)*nx + ip1 890 ! 891 src_lats(1) = src_center_lat(srch_add) 892 src_lats(2) = src_center_lat(e_add) 893 src_lats(3) = src_center_lat(ne_add) 894 src_lats(4) = src_center_lat(n_add) 895 ! 896 src_lons(1) = src_center_lon(srch_add) 897 src_lons(2) = src_center_lon(e_add) 898 src_lons(3) = src_center_lon(ne_add) 899 src_lons(4) = src_center_lon(n_add) 900 ! 901 !*** 902 !*** for consistency, we must make sure all lons are in 903 !*** same 2pi interval 904 !*** 905 ! 906 vec1_lon = src_lons(1) - plon 907 IF (vec1_lon > pi) THEN 908 src_lons(1) = src_lons(1) - pi2 909 ELSE IF (vec1_lon < -pi) THEN 910 src_lons(1) = src_lons(1) + pi2 911 ENDIF 912 DO n=2,4 913 vec1_lon = src_lons(n) - src_lons(1) 914 IF (vec1_lon > pi) THEN 915 src_lons(n) = src_lons(n) - pi2 916 ELSE IF (vec1_lon < -pi) THEN 917 src_lons(n) = src_lons(n) + pi2 918 ENDIF 919 END DO 920 ! 921 corner_loop: DO n=1,4 922 next_n = MOD(n,4) + 1 923 !*** 924 !*** here we take the cross product of the vector making 925 !*** up each box side with the vector formed by the vertex 926 !*** and search point. if all the cross products are 927 !*** positive, the point is contained in the box. 928 !*** 929 vec1_lat = src_lats(next_n) - src_lats(n) 930 vec1_lon = src_lons(next_n) - src_lons(n) 931 vec2_lat = plat - src_lats(n) 932 vec2_lon = plon - src_lons(n) 933 !*** 934 !*** check for 0,2pi crossings 935 !*** 936 IF (vec1_lon > three*pih) THEN 937 vec1_lon = vec1_lon - pi2 938 ELSE IF (vec1_lon < -three*pih) THEN 939 vec1_lon = vec1_lon + pi2 940 ENDIF 941 IF (vec2_lon > three*pih) THEN 942 vec2_lon = vec2_lon - pi2 943 ELSE IF (vec2_lon < -three*pih) THEN 944 vec2_lon = vec2_lon + pi2 945 ENDIF 946 ! 947 cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat 948 ! 949 !*** 950 !*** if cross product is less than zero, this cell 951 !*** doesn't work 952 !*** 953 IF (n == 1) cross_product_last = cross_product 954 IF (cross_product*cross_product_last < zero) & 955 EXIT corner_loop 956 cross_product_last = cross_product 957 ! 958 END DO corner_loop 959 !*** 960 !*** if cross products all same sign, we found the location 961 !*** 962 IF (n > 4) THEN 963 src_add(1) = srch_add 964 src_add(2) = e_add 965 src_add(3) = ne_add 966 src_add(4) = n_add 967 ! 968 lastsrc_add = srch_add 969 RETURN 970 ENDIF 971 !*** 972 !*** otherwise move on to next cell 973 !*** 974 ENDIF !bounding box check 975 END DO srch_loop 976 977 978 ENDDO 979 980 981 !*** 982 !*** if no cell found, point is likely either in a box that 983 !*** straddles either pole or is outside the grid. fall back 984 !*** to a distance-weighted average of the four closest 985 !*** points. go ahead and compute weights here, but store 986 !*** in src_lats and return -add to prevent the parent 987 !*** routine from computing bilinear weights 988 !*** 989 !print *,'Could not find location for ',plat,plon 990 !print *,'Using nearest-neighbor average for this point' 991 ! 992 coslat_dst = COS(plat) 993 sinlat_dst = SIN(plat) 994 coslon_dst = COS(plon) 995 sinlon_dst = SIN(plon) 996 ! 997 dist_min = bignum 998 src_lats = bignum 999 DO srch_add = min_add,max_add 1000 ! distance = ACOS(coslat_dst*COS(src_center_lat(srch_add))* & 1001 ! (coslon_dst*COS(src_center_lon(srch_add)) + & 1002 ! sinlon_dst*SIN(src_center_lon(srch_add)))+ & 1003 ! sinlat_dst*SIN(src_center_lat(srch_add))) 1004 1005 distance = ACOS(coslat_dst*cos_src_center_lat(srch_add)* & 1006 (coslon_dst*cos_src_center_lon(srch_add) + & 1007 sinlon_dst*sin_src_center_lon(srch_add))+ & 1008 sinlat_dst*sin_src_center_lat(srch_add)) 1009 1010 IF (distance < dist_min) THEN 1011 sort_loop: DO n=1,4 1012 IF (distance < src_lats(n)) THEN 1013 DO i=4,n+1,-1 1014 src_add (i) = src_add (i-1) 1015 src_lats(i) = src_lats(i-1) 1016 END DO 1017 src_add (n) = -srch_add 1018 src_lats(n) = distance 1019 dist_min = src_lats(4) 1020 EXIT sort_loop 1021 ENDIF 1022 END DO sort_loop 1023 ENDIF 1024 END DO 1025 ! 1026 src_lons = one/(src_lats + tiny) 1027 distance = SUM(src_lons) 1028 src_lats = src_lons/distance 1029 1030 !----------------------------------------------------------------------- 1031 1032 END SUBROUTINE grid_search_bilin 1033 1034 1035 !*********************************************************************** 1036 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1037 !************************************************************************ 1038 ! SUBROUTINE STORE_LINK_BILIN 1039 !************************************************************************ 1040 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1041 1042 SUBROUTINE store_link_bilin(dst_add, src_add, weights, nmap) 1043 1044 !----------------------------------------------------------------------- 1045 ! this routine stores the address and weight for four links 1046 ! associated with one destination point in the appropriate address 1047 ! and weight arrays and resizes those arrays if necessary. 1048 !----------------------------------------------------------------------- 1049 1050 !----------------------------------------------------------------------- 1051 ! input variables 1052 !----------------------------------------------------------------------- 1053 ! 1054 INTEGER :: dst_add,nmap 1055 ! 1056 INTEGER, DIMENSION(4) :: src_add 1057 ! 1058 REAL*8, DIMENSION(4) :: weights 1059 1060 !----------------------------------------------------------------------- 1061 ! 1062 ! local variables 1063 ! 1064 !----------------------------------------------------------------------- 1065 1066 INTEGER :: n,num_links_old 1067 1068 !----------------------------------------------------------------------- 1069 ! increment number of links and check to see if remap arrays need 1070 ! to be increased to accomodate the new link. then store the 1071 ! link. 1072 !----------------------------------------------------------------------- 1073 1074 num_links_old = num_links_map1 1075 num_links_map1 = num_links_old + 4 1076 1077 IF (num_links_map1 > max_links_map1) & 1078 CALL resize_remap_vars(1,resize_increment) 1079 1080 DO n=1,4 1081 grid1_add_map1(num_links_old+n) = src_add(n) 1082 grid2_add_map1(num_links_old+n) = dst_add 1083 wts_map1 (1,num_links_old+n) = weights(n) 1084 END DO 1085 1086 !----------------------------------------------------------------------- 1087 1088 END SUBROUTINE store_link_bilin 1089 1090 1091 1092 1093 1094 1095 1096 1097 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1098 !************************************************************************ 1099 ! SUBROUTINE INIT_REMAP_VARS 1100 !************************************************************************ 1101 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1102 ! 1103 SUBROUTINE init_remap_vars 1104 1105 !----------------------------------------------------------------------- 1106 ! 1107 ! this routine initializes some variables and provides an initial 1108 ! allocation of arrays (fairly large so frequent resizing 1109 ! unnecessary). 1110 ! 1111 !----------------------------------------------------------------------- 1112 1113 !----------------------------------------------------------------------- 1114 ! determine the number of weights 1115 !----------------------------------------------------------------------- 1116 num_wts = 1 ! bilinear interpolation 1117 !----------------------------------------------------------------------- 1118 ! initialize num_links and set max_links to four times the largest 1119 ! of the destination grid sizes initially (can be changed later). 1120 ! set a default resize increment to increase the size of link 1121 ! arrays if the number of links exceeds the initial size 1122 !----------------------------------------------------------------------- 1123 1124 num_links_map1 = 0 1125 max_links_map1 = 4*grid2_size 1126 IF (num_maps > 1) THEN 1127 num_links_map2 = 0 1128 max_links_map1 = MAX(4*grid1_size,4*grid2_size) 1129 max_links_map2 = max_links_map1 1130 ENDIF 1131 1132 resize_increment = 0.1*MAX(grid1_size,grid2_size) 1133 1134 !----------------------------------------------------------------------- 1135 ! allocate address and weight arrays for mapping 1 1136 !----------------------------------------------------------------------- 1137 1138 ALLOCATE (grid1_add_map1(max_links_map1), & 1139 grid2_add_map1(max_links_map1), & 1140 wts_map1(num_wts, max_links_map1)) 1141 1142 grid1_add_map1 = 0. 1143 grid2_add_map1 = 0. 1144 wts_map1 = 0. 1145 1146 !----------------------------------------------------------------------- 1147 1148 END SUBROUTINE init_remap_vars 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 !*********************************************************************** 1164 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1165 !************************************************************************ 1166 ! SUBROUTINE RESIZE_REMAP_VAR 1167 !************************************************************************ 1168 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1169 1170 SUBROUTINE resize_remap_vars(nmap, increment) 1171 1172 !----------------------------------------------------------------------- 1173 ! this routine resizes remapping arrays by increasing(decreasing) 1174 ! the max_links by increment 1175 !----------------------------------------------------------------------- 1176 1177 !----------------------------------------------------------------------- 1178 ! input variables 1179 !----------------------------------------------------------------------- 1180 1181 INTEGER :: & 1182 nmap, & ! identifies which mapping array to resize 1183 increment ! the number of links to add(subtract) to arrays 1184 1185 !----------------------------------------------------------------------- 1186 ! local variables 1187 !----------------------------------------------------------------------- 1188 1189 INTEGER :: & 1190 ierr, & ! error flag 1191 mxlinks ! size of link arrays 1192 1193 INTEGER, DIMENSION(:), ALLOCATABLE :: & 1194 add1_tmp, & ! temp array for resizing address arrays 1195 add2_tmp ! temp array for resizing address arrays 1196 ! 1197 ! temp array for resizing weight arrays 1198 ! 1199 REAL*8, DIMENSION(:,:), ALLOCATABLE :: wts_tmp 1200 ! 1201 !----------------------------------------------------------------------- 1202 !*** 1203 !*** allocate temporaries to hold original values 1204 !*** 1205 mxlinks = SIZE(grid1_add_map1) 1206 ALLOCATE (add1_tmp(mxlinks), add2_tmp(mxlinks), & 1207 wts_tmp(num_wts,mxlinks)) 1208 1209 add1_tmp = grid1_add_map1 1210 add2_tmp = grid2_add_map1 1211 wts_tmp = wts_map1 1212 1213 !*** 1214 !*** deallocate originals and increment max_links then 1215 !*** reallocate arrays at new size 1216 !*** 1217 1218 DEALLOCATE (grid1_add_map1, grid2_add_map1, wts_map1) 1219 max_links_map1 = mxlinks + increment 1220 ALLOCATE (grid1_add_map1(max_links_map1), & 1221 grid2_add_map1(max_links_map1), & 1222 wts_map1(num_wts,max_links_map1)) 1223 !*** 1224 !*** restore original values from temp arrays and 1225 !*** deallocate temps 1226 !*** 1227 mxlinks = MIN(mxlinks, max_links_map1) 1228 grid1_add_map1(1:mxlinks) = add1_tmp (1:mxlinks) 1229 grid2_add_map1(1:mxlinks) = add2_tmp (1:mxlinks) 1230 wts_map1 (:,1:mxlinks) = wts_tmp(:,1:mxlinks) 1231 DEALLOCATE(add1_tmp, add2_tmp, wts_tmp) 1232 1233 !----------------------------------------------------------------------- 1234 ! 1235 END SUBROUTINE resize_remap_vars 1236 ! 1237 !************************************************************************ 1238 ! 1239 END MODULE bilinear_interp 1240 1241 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1242 -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/dom_oce.F90
r12414 r13024 41 41 REAL(wp), PUBLIC :: e3zps_min !: miminum thickness for partial steps (meters) 42 42 REAL(wp), PUBLIC :: e3zps_rat !: minimum thickness ration for partial steps 43 INTEGER , PUBLIC :: nn_closea !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 43 44 44 45 INTEGER, PUBLIC :: nn_interp 45 CHARACTER(LEN= 132), PUBLIC :: cn_topo46 CHARACTER(LEN=200), PUBLIC :: cn_topo 46 47 CHARACTER(LEN=132), PUBLIC :: cn_bath 47 48 CHARACTER(LEN=132), PUBLIC :: cn_lon 48 49 CHARACTER(LEN=132), PUBLIC :: cn_lat 50 REAL(wp), PUBLIC :: rn_scale 49 51 50 52 LOGICAL, PUBLIC :: lzoom = .FALSE. !: zoom flag … … 279 281 #if defined key_agrif 280 282 LOGICAL, PUBLIC, PARAMETER :: lk_agrif = .TRUE. !: agrif flag 283 LOGICAL, PUBLIC :: lk_south, lk_north, lk_west, lk_east !: Child grid boundaries (interpolation or not) 281 284 #else 282 285 LOGICAL, PUBLIC, PARAMETER :: lk_agrif = .FALSE. !: agrif flag … … 326 329 & ff_f (jpi,jpj) , ff_t (jpi,jpj) , STAT=ierr(3) ) 327 330 ! 328 ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , STAT=ierr(4) ) 331 ! ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) , & 332 ! & gdept_b(jpi,jpj,jpk) , gdepw_b(jpi,jpj,jpk) , & 333 ! & gdept_n(jpi,jpj,jpk) , gdepw_n(jpi,jpj,jpk) , gde3w_n(jpi,jpj,jpk) , STAT=ierr(4) ) 334 335 ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk), STAT=ierr(4) ) 336 337 ! 338 ! ALLOCATE( e3t_0(jpi,jpj,jpk) , e3u_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , & 339 ! & e3t_b(jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b(jpi,jpj,jpk) , e3w_b(jpi,jpj,jpk) , & 340 ! & e3t_n(jpi,jpj,jpk) , e3u_n(jpi,jpj,jpk) , e3v_n(jpi,jpj,jpk) , e3f_n(jpi,jpj,jpk) , e3w_n(jpi,jpj,jpk) , & 341 ! & e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk) , & 342 ! ! ! 343 ! & e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , & 344 ! & e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) , & 345 ! & e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , STAT=ierr(5) ) 346 347 ! 348 ! ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , & 349 ! & hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) , & 350 ! & ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) , & 351 ! & hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6) ) 329 352 ! 330 353 ALLOCATE( e3t_0 (jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , & 331 & e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , STAT=ierr(5) ) 354 & e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , STAT=ierr(5) ) 332 355 ! 333 356 ALLOCATE( gdept_1d(jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) ,gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(6) ) 334 357 ! 335 ALLOCATE( bathy(jpi,jpj),mbathy(jpi,jpj), tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , & 358 ALLOCATE( bathy(jpi,jpj),mbathy(jpi,jpj), tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , & 336 359 & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , & 337 360 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(7) ) … … 340 363 & risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(8) ) 341 364 ! 342 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) , & 365 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) , & 343 366 & vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , wmask(jpi,jpj,jpk) , STAT=ierr(9) ) 344 367 ! … … 351 374 ALLOCATE( msk_opnsea (jpi,jpj), msk_csundef (jpi,jpj), & 352 375 & msk_csglo (jpi,jpj), msk_csrnf (jpi,jpj), msk_csemp (jpi,jpj), & 353 & msk_csgrpglo(jpi,jpj), msk_csgrprnf(jpi,jpj), msk_csgrpemp(jpi,jpj), STAT=ierr(11) ) 354 ! 376 & msk_csgrpglo(jpi,jpj), msk_csgrprnf(jpi,jpj), msk_csgrpemp(jpi,jpj), STAT=ierr(11) ) ! 355 377 dom_oce_alloc = MAXVAL(ierr) 356 378 ! -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domain.F90
r12870 r13024 62 62 !! - 1D configuration, move Coriolis, u and v at T-point 63 63 !!---------------------------------------------------------------------- 64 INTEGER :: jk ! dummy loop indices 65 INTEGER :: iconf = 0 ! local integers 66 REAL(wp), POINTER, DIMENSION(:,:) :: z1_hu_0, z1_hv_0 67 INTEGER , DIMENSION(jpi,jpj) :: ik_top , ik_bot ! top and bottom ocean level 68 !!---------------------------------------------------------------------- 69 ! 70 ! IF( nn_timing == 1 ) CALL timing_start('dom_init') 64 71 ! 65 72 IF(lwp) THEN … … 71 78 ! !== Reference coordinate system ==! 72 79 ! 73 CALL dom_nam ! read namelist ( namrun, namdom ) 74 ! 75 CALL dom_hgr ! Horizontal mesh 76 ! 77 CALL dom_zgr ! Vertical mesh and bathymetry 78 ! 79 CALL dom_msk ! compute mask (needed by write_cfg) 80 ! 81 IF ( ln_domclo ) CALL dom_clo ! Closed seas and lake 80 CALL dom_nam ! read namelist ( namrun, namdom ) 81 ! CALL dom_clo ! Closed seas and lake 82 83 CALL dom_hgr ! Horizontal mesh 84 CALL dom_zgr( ik_top, ik_bot ) ! Vertical mesh and bathymetry 85 CALL dom_msk( ik_top, ik_bot ) ! Masks 86 ! 82 87 ! 83 88 CALL dom_ctl ! print extrema of masked scale factors 84 89 ! 90 #if ! defined key_agrif 85 91 CALL cfg_write ! create the configuration file 92 #endif 86 93 ! 87 94 END SUBROUTINE dom_init … … 103 110 & nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, nn_euler , & 104 111 & ln_cfmeta, ln_iscpl 105 NAMELIST/namdom/ nn_bathy, cn_topo, cn_bath, cn_lon, cn_lat, nn_interp, & 106 & rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, & 107 & rn_atfp , rn_rdt , ln_crs , jphgr_msh , & 112 ! NAMELIST/namdom/ nn_bathy, cn_topo, cn_bath, cn_lon, cn_lat, nn_interp, & 113 ! & rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, & 114 ! & rn_atfp , rn_rdt , ln_crs , jphgr_msh , & 115 ! & ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, & 116 ! & ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 117 ! & ppa2, ppkth2, ppacr2 118 119 NAMELIST/namdom/ nn_bathy, cn_topo, cn_bath, cn_lon, cn_lat, rn_scale, nn_interp, & 120 & rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, & 121 & rn_atfp , rn_rdt , ln_crs , jphgr_msh , & 108 122 & ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, & 109 123 & ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 110 124 & ppa2, ppkth2, ppacr2 111 112 113 125 114 126 INTEGER :: ios ! Local integer output status for namelist read … … 181 193 182 194 183 195 ! 196 cn_topo ='' 197 cn_bath ='' 198 cn_lon ='' 199 cn_lat ='' 200 rn_scale = 1. 184 201 185 202 REWIND( numnam_ref ) ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep) … … 193 210 IF(lwm) WRITE ( numond, namdom ) 194 211 ! 212 213 214 195 215 IF(lwp) THEN 196 216 WRITE(numout,*) … … 198 218 WRITE(numout,*) ' flag read/compute bathymetry nn_bathy = ', nn_bathy 199 219 IF( nn_bathy == 2 ) THEN 200 WRITE(numout,*) ' compute bathymetry from file cn_topo = ', cn_topo 220 WRITE(numout,*) ' compute bathymetry from file cn_topo = ' , cn_topo 221 WRITE(numout,*) ' bathymetry name in file cn_bath = ' , cn_bath 222 WRITE(numout,*) ' longitude name in file cn_lon = ' , cn_lon 223 WRITE(numout,*) ' latitude name in file cn_lat = ' , cn_lat 224 WRITE(numout,*) ' bathmetry scale factor rn_scale = ' , rn_scale 201 225 ENDIF 202 226 WRITE(numout,*) ' Depth (if =0 bathy=jpkm1) rn_bathy = ', rn_bathy … … 255 279 !!---------------------------------------------------------------------- 256 280 ! 281 #undef CHECK_DOM 282 #ifdef CHECK_DOM 257 283 IF(lk_mpp) THEN 258 284 CALL mpp_minloc( 'dom_ctl', e1t(:,:), tmask_i(:,:), ze1min, iloc ) … … 292 318 WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 293 319 ENDIF 320 #endif 294 321 ! 295 322 ! check that all processes are still there... If some process have an error, -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/dombat.F90
r12414 r13024 2 2 3 3 USE dom_oce ! ocean domain 4 ! USE closea ! closed seas5 4 ! 6 5 USE in_out_manager ! I/O manager … … 8 7 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 9 8 USE lib_mpp ! distributed memory computing library 9 USE timing ! Timing 10 #if defined key_agrif 10 11 USE agrif_modutil 12 USE agrif_parameters 13 #endif 11 14 USE bilinear_interp 12 15 … … 20 23 SUBROUTINE dom_bat 21 24 22 INTEGER :: inum, isize, jsize, id, ji, jj 25 INTEGER :: inum, id, ji, jj,ji1,jj1 26 INTEGER :: iimin,iimax,jjmin,jjmax 23 27 INTEGER :: tabdim1, tabdim2, nxhr, nyhr, nxyhr 24 28 INTEGER, DIMENSION(2) :: ddims 25 29 INTEGER, DIMENSION(3) :: status 26 INTEGER, DIMENSION(:,:), ALLOCATABLE :: trouble_points , vardep,mask_oce27 INTEGER :: iimin,iimax,jjmin,jjmax28 30 INTEGER, DIMENSION(1) :: i_min,i_max 29 INTEGER, DIMENSION(1) ::j_min,j_max 30 REAL(wp), DIMENSION(jpi,jpj) :: myglamf 31 INTEGER,DIMENSION(:) ,POINTER :: src_add,dst_add => NULL() 32 REAL(wp), DIMENSION(:) ,ALLOCATABLE :: vardep1d, lon_new1D,lat_new1D 33 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bathy_new, lat_new, lon_new, bathy_test 34 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: coarselon, coarselat, coarsebathy 35 REAL(wp) :: Cell_lonmin, Cell_lonmax, Cell_latmin, Cell_latmax 31 INTEGER, DIMENSION(1) :: j_min,j_max 32 INTEGER, DIMENSION(:) ,ALLOCATABLE :: src_add,dst_add 33 INTEGER, DIMENSION(:,:), ALLOCATABLE :: trouble_points , vardep,mask_oce 34 REAL(wp) ,DIMENSION(jpi,jpj):: Cell_lonmin, Cell_lonmax, Cell_latmin, Cell_latmax 35 REAL(wp) ::zdel 36 REAL(wp), DIMENSION(:) , ALLOCATABLE :: lon_new1D , lat_new1D, vardep1d 37 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: coarselon, coarselat, coarsebathy, bathy_test 38 REAL(wp), DIMENSION(:,:),ALLOCATABLE :: matrix,interpdata 39 LOGICAL :: lonlat_2D, ln_pacifiq 36 40 LOGICAL :: identical_grids 37 41 LOGICAL, DIMENSION(:,:), ALLOCATABLE :: masksrc 38 REAL *8, DIMENSION(:,:),POINTER :: matrix,interpdata => NULL()39 LOGICAL :: lonlat_2D40 42 REAL(wp), DIMENSION(jpi,jpj) :: zglamt, zgphi, zglamu, zglamv, zglamf 43 REAL(wp) :: zshift 44 41 45 CHARACTER(32) :: bathyfile, bathyname, lonname,latname 42 43 lonlat_2D=.false.44 46 45 47 bathyfile=TRIM(cn_topo) … … 49 51 50 52 CALL iom_open( bathyfile, inum, lagrif=.FALSE. ) 53 54 ! check if lon/lat are 2D arrays 55 id = iom_varid( inum, lonname, ddims ) 56 IF (ddims(2)==0) THEN 57 lonlat_2D = .FALSE. 58 ELSE 59 lonlat_2D = .TRUE. 60 ENDIF 61 51 62 id = iom_varid( inum, bathyname, ddims ) 52 63 ln_pacifiq = .FALSE. 64 zglamt(:,:) = glamt(:,:) 65 zglamu(:,:) = glamu(:,:) 66 zglamv(:,:) = glamv(:,:) 67 zglamf(:,:) = glamf(:,:) 68 69 IF( glamt(1,1) .GT. glamt(jpi,jpj) ) ln_pacifiq =.false. 70 71 zshift = 0. 72 IF( ln_pacifiq ) THEN 73 zshift = 0.!Abs(minval(glamt)) +0.1 74 WHERE ( glamt < 0 ) 75 zglamt = zglamt + zshift + 360. 76 END WHERE 77 WHERE ( glamu < 0 ) 78 zglamu = zglamu + zshift +360. 79 END WHERE 80 WHERE ( glamv < 0 ) 81 zglamv = zglamv + zshift +360. 82 END WHERE 83 WHERE ( glamf < 0 ) 84 zglamf = zglamf + zshift +360. 85 END WHERE 86 ENDIF 87 53 88 status=-1 54 ALLOCATE(lon_new (ddims(1),ddims(2)), STAT=status(1))55 ALLOCATE(lat_new (ddims(1),ddims(2)), STAT=status(2))56 ALLOCATE(bathy_new(ddims(1),ddims(2)), STAT=status(3))57 IF( sum(status) /= 0 ) CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' )58 89 59 90 IF (lonlat_2D) THEN 60 CALL iom_get ( inum, jpdom_unknown, lonname, lon_new ) 61 CALL iom_get ( inum, jpdom_unknown, latname, lat_new ) 91 ! here implicitly it's old topo database (orca format) 92 ALLOCATE(coarselon (ddims(1),ddims(2)), STAT=status(1)) 93 ALLOCATE(coarselat (ddims(1),ddims(2)), STAT=status(2)) 94 ALLOCATE(coarsebathy(ddims(1),ddims(2)), STAT=status(3)) 95 IF( sum(status) /= 0 ) CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' ) 96 CALL iom_get ( inum, jpdom_unknown, lonname, coarselon ) 97 CALL iom_get ( inum, jpdom_unknown, latname, coarselat ) 98 CALL iom_get ( inum, jpdom_unknown, bathyname, coarsebathy ) 99 CALL iom_close (inum) 100 IF( ln_pacifiq ) THEN 101 ! WHERE(coarselon < 0.00001) 102 coarselon = Coarselon + zshift 103 ! END WHERE 104 ENDIF 105 ! equivalent to new database 62 106 ELSE 63 ALLOCATE(lon_new1D(ddims(1)), lat_new1D(ddims(2))) 64 CALL iom_get ( inum, jpdom_unknown, lonname, lon_new1D ) 65 CALL iom_get ( inum, jpdom_unknown, latname, lat_new1D ) 66 DO ji=1, ddims(1) 67 lon_new(ji,:)=lon_new1D(ji) 68 ENDDO 69 DO ji=1, ddims(2) 70 lat_new(:,ji)=lat_new1D(ji) 71 ENDDO 107 ALLOCATE(lon_new1D(ddims(1)), lat_new1D(ddims(2))) 108 CALL iom_get ( inum, jpdom_unknown, lonname, lon_new1D ) 109 CALL iom_get ( inum, jpdom_unknown, latname, lat_new1D ) 110 IF( ln_pacifiq ) THEN 111 WHERE(lon_new1D < 0.00001) 112 lon_new1D = lon_new1D +360.!zshift 113 END WHERE 114 ENDIF 115 zdel = 0.00 116 IF( MAXVAL(zglamf) > 180 + zshift ) THEN 117 ! 118 ! WHERE( lon_new1D < 0 ) 119 ! lon_new1D = lon_new1D + 360. 120 ! END WHERE 121 ! 122 i_min = MAXLOC(lon_new1D,mask = lon_new1D < MINVAL(zglamf(1:jpi-1,1:jpj-1)) ) 123 i_max = MINLOC(lon_new1D,mask = lon_new1D > MAXVAL(zglamf(1:jpi-1,1:jpj-1)) ) 124 j_min = MAXLOC(lat_new1D,mask = lat_new1D < MINVAL( gphif(1:jpi-1,1:jpj-1)) ) 125 j_max = MINLOC(lat_new1D,mask = lat_new1D > MAXVAL( gphif(1:jpi-1,1:jpj-1)) ) 126 ! 127 tabdim1 = ( SIZE(lon_new1D) - i_min(1) + 1 ) + i_max(1) 128 ! 129 IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(lat_new1D,1) ) THEN 130 j_min(1) = j_min(1) - 2 131 j_max(1) = j_max(1)+ 3 132 ENDIF 133 tabdim2 = j_max(1) - j_min(1) + 1 134 ! 135 ALLOCATE(coarselon (tabdim1,tabdim2), STAT=status(1)) 136 ALLOCATE(coarselat (tabdim1,tabdim2), STAT=status(2)) 137 ALLOCATE(Coarsebathy(tabdim1,tabdim2), STAT=status(3)) 138 139 IF( SUM(status) /= 0 ) CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' ) 140 ! 141 DO ji = 1,tabdim1 142 coarselat(ji,:) = lat_new1D(j_min(1):j_max(1)) 143 END DO 144 145 ! 146 DO jj = 1, tabdim2 147 coarselon(1:SIZE(lon_new1D)-i_min(1)+1 ,jj) = lon_new1D(i_min(1):SIZE(lon_new1D)) 148 coarselon(2+SIZE(lon_new1D)-i_min(1):tabdim1,jj) = lon_new1D(1:i_max(1)) 149 END DO 150 ! 151 CALL iom_get(inum, jpdom_unknown, bathyname,coarsebathy(1:SIZE(lon_new1D)-i_min(1)+1,:), & 152 kstart=(/i_min(1),j_min(1)/), kcount=(/SIZE(lon_new1D)-i_min(1)+1,tabdim2/)) ! +1? 153 ! 154 CALL iom_get(inum, jpdom_unknown, bathyname,coarsebathy(2+SIZE(lon_new1D)-i_min(1):tabdim1,:), & 155 kstart=(/1,j_min(1)/),kcount=(/i_max(1),tabdim2/)) 156 ! 157 ELSE 158 ! WHERE( lon_new1D > (180. + zshift) ) lon_new1D = lon_new1D - 360. 159 i_min = MAXLOC(lon_new1D,mask = lon_new1D < MINVAL(zglamf)-zdel) 160 i_max = MINLOC(lon_new1D,mask = lon_new1D > MAXVAL(zglamf)+zdel) 161 j_min = MAXLOC(lat_new1D,mask = lat_new1D < MINVAL(gphif)-zdel) 162 j_max = MINLOC(lat_new1D,mask = lat_new1D > MAXVAL(gphif)+zdel) 163 164 i_min(1)=max(i_min(1),1) 165 ! 166 IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(lon_new1D,1) ) THEN 167 i_min(1) = i_min(1)-2 168 i_max(1) = i_max(1)+3 169 ENDIF 170 tabdim1 = i_max(1) - i_min(1) + 1 171 ! 172 IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(lat_new1D,1) ) THEN 173 j_min(1) = j_min(1)-2 174 j_max(1) = j_max(1)+3 175 ENDIF 176 tabdim2 = j_max(1) - j_min(1) + 1 177 ! 178 ALLOCATE(coarselon (tabdim1,tabdim2), STAT=status(1)) 179 ALLOCATE(coarselat (tabdim1,tabdim2), STAT=status(2)) 180 ALLOCATE(coarsebathy(tabdim1,tabdim2), STAT=status(3)) 181 182 IF( SUM(status) /= 0 ) CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' ) 183 ! 184 DO jj = 1,tabdim2 185 coarselon(:,jj) = lon_new1D(i_min(1):i_max(1)) 186 END DO 187 ! 188 DO ji = 1,tabdim1 189 coarselat(ji,:) = lat_new1D(j_min(1):j_max(1)) 190 END DO 191 ! 192 CALL iom_get(inum,jpdom_unknown,bathyname,coarsebathy, & 193 & kstart=(/i_min(1),j_min(1)/),kcount=(/tabdim1,tabdim2/)) 194 ! 195 ENDIF ! > 180 196 197 DEALLOCATE(lon_new1D) ; DEALLOCATE(lat_new1D) 198 CALL iom_close (inum) 199 200 coarsebathy = coarsebathy *rn_scale 201 ! reset land to 0 202 WHERE (coarsebathy < 0.) 203 coarsebathy=0. 204 ENDWHERE 205 206 ENDIF ! external 207 208 IF(lwp) THEN 209 WRITE(numout,*) 'Interpolation of high resolution bathymetry on child grid' 210 IF( nn_interp == 0 ) THEN 211 WRITE(numout,*) 'Arithmetic average ...' 212 ELSE IF( nn_interp == 1 ) THEN 213 WRITE(numout,*) 'Median average ...' 214 ELSE IF( nn_interp == 2 ) THEN 215 WRITE(numout,*) 'Bilinear interpolation ...' 216 ELSE 217 WRITE(*,*) 'bad value for nn_interp variable ( must be 0,1 or 2 )' 218 STOP 219 ENDIF 72 220 ENDIF 73 CALL iom_get ( inum, jpdom_unknown, bathyname, bathy_new ) 74 CALL iom_close (inum) 75 WHERE (bathy_new > 0.) 76 bathy_new=0. 77 ENDWHERE 78 bathy_new=-bathy_new 79 80 ! Eventually add here a pre-selection of the area (coarselon/lat) 81 i_min=10 ; j_min=10 82 i_max= ddims(1)-10 ; j_max=ddims(2)-10 83 84 tabdim1 = i_max(1) - i_min(1) + 1 85 tabdim2 = j_max(1) - j_min(1) + 1 86 ! 87 88 ALLOCATE(coarselon(tabdim1,tabdim2)) 89 ALLOCATE(coarselat(tabdim1,tabdim2)) 90 ALLOCATE(coarsebathy(tabdim1,tabdim2)) 91 ! 92 WHERE( lon_new < 0. ) lon_new = lon_new + 360. 93 myglamf=glamf 94 WHERE( myglamf < 0. ) myglamf = myglamf + 360. 95 96 coarselat(:,:) = lat_new (i_min(1):i_max(1), j_min(1):j_max(1)) 97 coarselon (:,:) = lon_new (i_min(1):i_max(1), j_min(1):j_max(1)) 98 coarsebathy(:,:) = bathy_new(i_min(1):i_max(1), j_min(1):j_max(1)) 99 100 IF( nn_interp == 0 .OR. nn_interp == 1 ) THEN ! arithmetic or median averages 101 ! ! ----------------------------- 102 ALLOCATE(trouble_points(jpi,jpj)) 103 trouble_points(:,:) = 0 104 ! 105 DO jj = 2, jpj 106 DO ji = 2, jpi 107 ! 108 ! fine grid cell extension 109 Cell_lonmin = MIN(myglamf(ji-1,jj-1),myglamf(ji,jj-1),myglamf(ji,jj),myglamf(ji-1,jj)) 110 Cell_lonmax = MAX(myglamf(ji-1,jj-1),myglamf(ji,jj-1),myglamf(ji,jj),myglamf(ji-1,jj)) 111 Cell_latmin = MIN(gphif(ji-1,jj-1),gphif(ji,jj-1),gphif(ji,jj),gphif(ji-1,jj)) 112 Cell_latmax = MAX(gphif(ji-1,jj-1),gphif(ji,jj-1),gphif(ji,jj),gphif(ji-1,jj)) 113 ! 114 ! look for points in G0 (bathy dataset) contained in the fine grid cells 115 iimin = 1 116 DO WHILE( coarselon(iimin,1) < Cell_lonmin ) 117 iimin = iimin + 1 118 ENDDO 119 ! 120 jjmin = 1 121 DO WHILE( coarselat(iimin,jjmin) < Cell_latmin ) 122 jjmin = jjmin + 1 123 ENDDO 124 ! 125 iimax = iimin 126 DO WHILE( coarselon(iimax,1) <= Cell_lonmax ) 127 iimax = iimax + 1 128 iimax = MIN( iimax,SIZE(coarsebathy,1)) 129 ENDDO 130 ! 131 jjmax = jjmin 132 DO WHILE( coarselat(iimax,jjmax) <= Cell_latmax ) 133 jjmax = jjmax + 1 134 jjmax = MIN( jjmax,SIZE(coarsebathy,2)) 135 ENDDO 136 ! 137 IF( .NOT. Agrif_Root() ) THEN 138 iimax = iimax-1 139 jjmax = jjmax-1 140 ELSE 141 iimax = MAX(iimin,iimax-1) 142 jjmax = MAX(jjmin,jjmax-1) 143 ENDIF 144 ! 145 iimin = MAX( iimin,1 ) 146 jjmin = MAX( jjmin,1 ) 147 iimax = MIN( iimax,SIZE(coarsebathy,1)) 148 jjmax = MIN( jjmax,SIZE(coarsebathy,2)) 149 150 nxhr = iimax - iimin + 1 151 nyhr = jjmax - jjmin + 1 152 153 154 IF( nxhr == 0 .OR. nyhr == 0 ) THEN 155 ! 156 trouble_points(ji,jj) = 1 157 ! 158 ELSE 159 ! 160 ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) ) 161 vardep(:,:) = coarsebathy(iimin:iimax,jjmin:jjmax) 162 ! 163 WHERE( vardep(:,:) .GT. 0. ) ; mask_oce = 1 ; 164 ELSEWHERE ; mask_oce = 0 ; 165 ENDWHERE 166 ! 167 nxyhr = nxhr*nyhr 168 IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0 169 bathy(ji,jj) = 0. 170 ELSE 171 IF( nn_interp == 0 ) THEN ! Arithmetic average 172 bathy(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) ) 173 ELSE ! Median average 174 ALLOCATE(vardep1d(nxyhr)) 175 vardep1d = RESHAPE(vardep,(/ nxyhr /) ) 176 !!CALL ssort(vardep1d,nxyhr) 177 CALL quicksort(vardep1d,1,nxyhr) 178 ! 179 ! Calculate median 180 IF (MOD(nxyhr,2) .NE. 0) THEN 181 bathy(ji,jj) = vardep1d( nxyhr/2 + 1 ) 182 ELSE 183 bathy(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) ) 184 END IF 185 DEALLOCATE(vardep1d) 186 ENDIF 187 ENDIF 188 DEALLOCATE (mask_oce,vardep) 189 ! 190 ENDIF 191 ENDDO 192 ENDDO 193 194 IF( SUM( trouble_points ) > 0 ) THEN 195 PRINT*,'too much empty cells, proceed to bilinear interpolation' 196 nn_interp = 2 197 stop 198 ENDIF 199 ENDIF 200 201 #undef MYTEST 202 #ifdef MYTEST 203 IF( nn_interp == 2) THEN ! Bilinear interpolation 204 ! ! ----------------------------- 205 identical_grids = .FALSE. 206 207 IF( SIZE(coarselat,1) == jpi .AND. SIZE(coarselat,2) == jpj .AND. & 208 SIZE(coarselon,1) == jpj .AND. SIZE(coarselon,2) == jpj ) THEN 209 IF( MAXVAL( ABS(coarselat(:,:)- gphit(:,:)) ) < 0.0001 .AND. & 210 MAXVAL( ABS(coarselon(:,:)- glamt(:,:)) ) < 0.0001 ) THEN 211 PRINT*,'same grid between ', cn_topo,' and child domain' 212 bathy = bathy_new 213 identical_grids = .TRUE. 214 ENDIF 215 ENDIF 216 217 IF( .NOT. identical_grids ) THEN 218 219 ALLOCATE(masksrc(SIZE(bathy_new,1),SIZE(bathy_new,2))) 220 ALLOCATE(bathy_test(jpi,jpj)) 221 ! 222 !Where(G0%bathy_meter.le.0.00001) 223 ! masksrc = .false. 224 !ElseWhere 225 masksrc = .TRUE. 226 !End where 227 ! 228 ! compute remapping matrix thanks to SCRIP package 229 CALL get_remap_matrix(coarselat,gphit,coarselon,glamt,masksrc,matrix,src_add,dst_add) 230 CALL make_remap(bathy_new,bathy_test,jpi,jpj,matrix,src_add,dst_add) 231 ! 232 bathy = bathy_test 233 ! 234 DEALLOCATE(masksrc) 235 DEALLOCATE(bathy_test) 236 237 ENDIF 221 bathy(:,:) = 0. 222 ! 223 !------------------------------------ 224 !MEDIAN AVERAGE or ARITHMETIC AVERAGE 225 !------------------------------------ 226 ! 227 IF( nn_interp == 0 .OR. nn_interp == 1 ) THEN 228 ! 229 ALLOCATE(trouble_points(jpi,jpj)) 230 trouble_points = 0 231 ! 232 ! POINT DETECTION 233 ! 234 ! 235 DO jj = 1,jpj 236 DO ji = 1,jpi 237 ! 238 ! FINE GRID CELLS DEFINITION 239 ji1=max(ji-1,1) 240 jj1=max(jj-1,1) 241 242 ! 243 Cell_lonmin(ji,jj) = MIN(zglamf(ji1,jj1),zglamf(ji,jj1),zglamf(ji,jj),zglamf(ji1,jj)) 244 Cell_lonmax(ji,jj) = MAX(zglamf(ji1,jj1),zglamf(ji,jj1),zglamf(ji,jj),zglamf(ji1,jj)) 245 Cell_latmin(ji,jj) = MIN(gphif(ji1,jj1),gphif(ji,jj1),gphif(ji,jj),gphif(ji1,jj)) 246 Cell_latmax(ji,jj) = MAX(gphif(ji1,jj1),gphif(ji,jj1),gphif(ji,jj),gphif(ji1,jj)) 247 IF( ABS(Cell_lonmax(ji,jj) - Cell_lonmin(ji,jj) ) > 180 ) THEN 248 zdel = Cell_lonmax(ji,jj) 249 Cell_lonmax(ji,jj) = Cell_lonmin(ji,jj) 250 Cell_lonmin(ji,jj) = zdel-360 251 ENDIF 252 ! 253 ! SEARCH FOR ALL POINTS CONTAINED IN THIS CELL 254 ! 255 ! ENDDO 256 ! ENDDO 257 ! CALL lbc_lnk( 'dom_bat', Cell_lonmin, 'T', 1. ) 258 ! CALL lbc_lnk( 'dom_bat', Cell_lonmax, 'T', 1. ) 259 ! CALL lbc_lnk( 'dom_bat', Cell_latmin, 'T', 1. ) 260 ! CALL lbc_lnk( 'dom_bat', Cell_latmax, 'T', 1. ) 261 262 263 ! DO jj = 2,jpj 264 ! DO ji = 2,jpi 265 iimin = 1 266 DO WHILE( coarselon(iimin,1) < Cell_lonmin(ji,jj) ) 267 iimin = iimin + 1 268 ! IF ( iimin .LE. 1 ) THEN 269 ! iimin = 1 270 ! EXIT 271 ! ENDIF 272 ENDDO 273 ! 274 jjmin = 1 275 DO WHILE( coarselat(iimin,jjmin) < Cell_latmin(ji,jj) ) 276 jjmin = jjmin + 1 277 ! IF ( iimin .LE. 1 ) THEN 278 ! iimin = 1 279 ! EXIT 280 ! ENDIF 281 ENDDO 282 jjmin=max(1,jjmin) 283 ! 284 iimax = iimin 285 DO WHILE( coarselon(iimax,1)<= Cell_lonmax(ji,jj) ) 286 iimax = iimax + 1 287 IF ( iimax .GE. SIZE(coarsebathy,1) ) THEN 288 iimax = MIN( iimax,SIZE(coarsebathy,1)) 289 EXIT 290 ENDIF 291 ENDDO 292 ! 293 jjmax = jjmin 294 DO WHILE( coarselat(iimax,jjmax) <= Cell_latmax(ji,jj) ) 295 jjmax = jjmax + 1 296 IF ( jjmax .GE. SIZE(coarsebathy,2) ) THEN 297 jjmax = MIN( jjmax,SIZE(coarsebathy,2)) 298 EXIT 299 ENDIF 300 ENDDO 301 ! 302 ! iimax = iimax-1 303 ! jjmax = jjmax-1 304 ! 305 iimin = MAX( iimin,1 ) 306 jjmin = MAX( jjmin,1 ) 307 iimax = MIN( iimax,SIZE(coarsebathy,1)) 308 jjmax = MIN( jjmax,SIZE(coarsebathy,2)) 309 310 nxhr = iimax - iimin + 1 311 nyhr = jjmax - jjmin + 1 312 313 314 315 IF( nxhr == 0 .OR. nyhr == 0 ) THEN 316 trouble_points(ji,jj) = 1 317 ELSE 318 ALLOCATE( vardep(nxhr,nyhr) ) 319 ALLOCATE( mask_oce(nxhr,nyhr) ) 320 mask_oce = 0 321 vardep(:,:) = coarsebathy(iimin:iimax,jjmin:jjmax) 322 WHERE( vardep(:,:) .GT. 0. ) 323 mask_oce = 1 324 ENDWHERE 325 nxyhr = nxhr*nyhr 326 ! IF( SUM(mask_oce) == 0 ) THEN 327 IF( SUM(mask_oce) < 0.5*(nxhr*nyhr) ) THEN 328 bathy(ji,jj) = 0. 329 ELSE 330 IF( nn_interp == 0 ) THEN 331 ! Arithmetic average 332 bathy(ji,jj) = SUM (vardep(:,:)*mask_oce(:,:))/SUM(mask_oce) 333 ELSE 334 ! Median average 335 ! 336 ALLOCATE(vardep1d(nxhr*nyhr)) 337 vardep1d = RESHAPE(vardep,(/ nxhr*nyhr /) ) 338 ! CALL ssort(vardep1d,nxhr*nyhr) 339 CALL quicksort(vardep1d,1,nxhr*nyhr) 340 ! 341 ! Calculate median 342 ! 343 IF (MOD(SUM(mask_oce),2) .NE. 0) THEN 344 bathy(ji,jj) = vardep1d( nxyhr/2 + 1 ) 345 ELSE 346 bathy(ji,jj) =0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) ) 347 END IF 348 DEALLOCATE(vardep1d) 349 ENDIF 350 ENDIF 351 DEALLOCATE (mask_oce,vardep) 352 ENDIF 353 ENDDO 354 ENDDO 355 356 IF( SUM( trouble_points ) > 0 ) THEN 357 CALL ctl_warn ('too much empty cells, proceed to bilinear interpolation') 358 nn_interp = 2 359 ENDIF 360 ENDIF 361 362 ! 363 ! create logical array masksrc 364 ! 365 IF( nn_interp == 2) THEN 366 ! 367 identical_grids = .FALSE. 368 369 # ifdef key_agrif 370 IF( Agrif_Parent(jpi) == jpi .AND. & 371 Agrif_Parent(jpj) == jpj ) THEN 372 IF( MAXVAL( ABS(coarselat(:,:)- gphit(:,:)) ) < 0.0001 .AND. & 373 MAXVAL( ABS(coarselon(:,:)- gphit(:,:)) ) < 0.0001 ) THEN 374 bathy(:,:) = coarsebathy(:,:) 375 IF(lwp) WRITE(numout,*) 'same grid between ', bathyname,' and child domain' 376 identical_grids = .TRUE. 377 ENDIF 378 ENDIF 379 # endif 380 IF( .NOT. identical_grids ) THEN 381 ALLOCATE(masksrc(SIZE(coarsebathy,1),SIZE(coarsebathy,2))) 382 ALLOCATE(bathy_test(jpi,jpj)) 383 ! 384 ! Where(G0%bathy_meter.le.0.00001) 385 ! masksrc = .false. 386 ! ElseWhere 387 ! 388 masksrc = .TRUE. 389 ! 390 ! End where 391 ! 392 ! compute remapping matrix thanks to SCRIP package 393 ! 394 CALL get_remap_matrix(coarselat,gphit, & 395 coarselon,zglamt, & 396 masksrc,matrix,src_add,dst_add) 397 CALL make_remap(coarsebathy,bathy_test,jpi,jpj, & 398 matrix,src_add,dst_add) 399 ! 400 bathy= bathy_test 401 ! 402 DEALLOCATE(masksrc) 403 DEALLOCATE(bathy_test) 404 ENDIF 238 405 ! 239 ENDIF 240 #endif 406 ENDIF 407 CALL lbc_lnk( 'dom_bat', bathy, 'T', 1. ) 408 409 ! Correct South and North 410 #if defined key_agrif 411 IF( ln_bry_south ) THEN 412 IF( (nbondj == -1).OR.(nbondj == 2) ) THEN 413 bathy(:,1)=bathy(:,2) 414 ENDIF 415 ELSE 416 bathy(:,1) = 0. 417 ENDIF 418 #else 419 IF ((nbondj == -1).OR.(nbondj == 2)) THEN 420 bathy(:,1)=bathy(:,2) 421 ENDIF 422 #endif 423 424 IF ((nbondj == 1).OR.(nbondj == 2)) THEN 425 bathy(:,jpj)=bathy(:,jpj-1) 426 ENDIF 427 428 ! Correct West and East 429 IF (jperio /=1) THEN 430 IF ((nbondi == -1).OR.(nbondi == 2)) THEN 431 bathy(1,:)=bathy(2,:) 432 ENDIF 433 IF ((nbondi == 1).OR.(nbondi == 2)) THEN 434 bathy(jpi,:)=bathy(jpi-1,:) 435 ENDIF 436 ENDIF 437 438 241 439 END SUBROUTINE dom_bat 242 440 -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domclo.F90
r12442 r13024 82 82 LOGICAL :: lskip ! flag in case lake seed on land or already filled (...) 83 83 84 NAMELIST/namclo/ rn_lon_opnsea, rn_lat_opnsea, nn_closea, sn_lake 85 84 86 !!---------------------------------------------------------------------- 85 87 !! 0 : Read namelist for closed sea definition … … 89 91 IF(lwp) WRITE(numout,*)'~~~~~~~' 90 92 91 NAMELIST/namclo/ rn_lon_opnsea, rn_lat_opnsea, nn_closea, sn_lake92 93 !!--------------------------------------------------------------------- 93 94 -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domhgr.F90
r12414 r13024 434 434 ! ------------------------------------------ 435 435 ! The equator line must be the latitude coordinate axe 436 !(PM) be carefull with nperio/jperio436 !(PM) be carefull with nperio/jperio 437 437 IF( jperio == 2 ) THEN 438 438 znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) … … 455 455 ! 456 456 INTEGER :: inum ! temporary logical unit 457 CHARACTER(LEN=135) :: coordinate_filename 457 458 !!---------------------------------------------------------------------- 458 459 ! … … 463 464 ENDIF 464 465 ! 465 CALL iom_open( 'coordinates', inum ) 466 IF (ln_read_cfg) THEN 467 coordinate_filename=TRIM(cn_domcfg) 468 ELSE 469 coordinate_filename='coordinates' 470 ENDIF 471 CALL iom_open( coordinate_filename, inum ) 466 472 ! 467 473 CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/dommsk.F90
r12414 r13024 26 26 USE domisf ! domain: ice shelf 27 27 USE domwri ! domain: write the meshmask file 28 USE usrdef_fmask ! user defined fmask 28 29 USE bdy_oce ! open boundary 29 30 ! … … 52 53 CONTAINS 53 54 54 SUBROUTINE dom_msk 55 SUBROUTINE dom_msk( k_top, k_bot ) 55 56 !!--------------------------------------------------------------------- 56 57 !! *** ROUTINE dom_msk *** … … 62 63 !! and ko_bot, the indices of the fist and last ocean t-levels which 63 64 !! are either defined in usrdef_zgr or read in zgr_read. 64 !! The velocity masks (umask, vmask )65 !! The velocity masks (umask, vmask, wmask, wumask, wvmask) 65 66 !! are deduced from a product of the two neighboring tmask. 66 67 !! The vorticity mask (fmask) is deduced from tmask taking … … 77 78 !! due to cyclic or North Fold boundaries as well as MPP halos. 78 79 !! 79 !! ** Action : tmask, umask, vmask, wmask : land/ocean mask80 !! ** Action : tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask 80 81 !! at t-, u-, v- w, wu-, and wv-points (=0. or 1.) 81 82 !! fmask : land/ocean mask at f-point (=0., or =1., or … … 85 86 !! ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask 86 87 !!---------------------------------------------------------------------- 88 89 INTEGER, DIMENSION(:,:), INTENT(in) :: k_top, k_bot ! first and last ocean level 87 90 ! 88 91 INTEGER :: ji, jj, jk ! dummy loop indices … … 133 136 ! N.B. tmask has already the right boundary conditions since mbathy is ok 134 137 ! 138 ! tmask(:,:,:) = 0._wp 139 ! DO jk = 1, jpk 140 ! DO jj = 1, jpj 141 ! DO ji = 1, jpi 142 ! IF( ( REAL( mbathy (ji,jj) - jk, wp ) + 0.1_wp >= 0._wp ) & 143 ! & .AND. ( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp <= 0._wp ) ) THEN 144 ! tmask(ji,jj,jk) = 1._wp 145 ! END IF 146 ! END DO 147 ! END DO 148 ! END DO 149 150 ! IF ( ln_isfsubgl ) CALL zgr_isf_subgl 151 152 ! Ocean/land mask at t-point (computed from ko_top and ko_bot) 153 ! ---------------------------- 154 ! 135 155 tmask(:,:,:) = 0._wp 136 DO jk = 1, jpk156 IF( ln_read_cfg) THEN 137 157 DO jj = 1, jpj 138 158 DO ji = 1, jpi 139 IF( ( REAL( mbathy (ji,jj) - jk, wp ) + 0.1_wp >= 0._wp ) & 140 & .AND. ( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp <= 0._wp ) ) THEN 141 tmask(ji,jj,jk) = 1._wp 142 END IF 159 iktop = k_top(ji,jj) 160 ikbot = k_bot(ji,jj) 161 IF( iktop /= 0 ) THEN ! water in the column 162 tmask(ji,jj,iktop:ikbot ) = 1._wp 163 ENDIF 164 END DO 165 END DO 166 ELSE 167 DO jk = 1, jpk 168 DO jj = 1, jpj 169 DO ji = 1, jpi 170 IF( ( REAL( mbathy (ji,jj) - jk, wp ) + 0.1_wp >= 0._wp ) & 171 & .AND. ( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp <= 0._wp ) ) THEN 172 tmask(ji,jj,jk) = 1._wp 173 END IF 174 END DO 143 175 END DO 144 176 END DO 145 END DO146 147 IF ( ln_isfsubgl ) CALL zgr_isf_subgl 177 IF ( ln_isfsubgl ) CALL zgr_isf_subgl 178 ENDIF 179 148 180 149 181 !SF add here lbc_lnk: bug not still understood : cause now domain configuration is read ! … … 272 304 #if defined key_agrif 273 305 IF( .NOT. AGRIF_Root() ) THEN 274 IF ((nbondi == 1).OR.(nbondi == 2))fmask(nlci-1 , : ,jk) = 0.e0 ! east275 IF ((nbondi == -1).OR.(nbondi == 2))fmask(1 , : ,jk) = 0.e0 ! west276 IF ((nbondj == 1).OR.(nbondj == 2)) fmask(: ,nlcj-1 ,jk) = 0.e0 ! north277 IF ((nbondj == -1).OR.(nbondj == 2)) fmask(: ,1 ,jk) = 0.e0 ! south306 IF(lk_east) fmask(nlci-1 , : ,jk) = 0.e0 ! east 307 IF(lk_west) fmask(1 , : ,jk) = 0.e0 ! west 308 IF(lk_north) fmask(: ,nlcj-1 ,jk) = 0.e0 ! north 309 IF(lk_south) fmask(: ,1 ,jk) = 0.e0 ! south 278 310 ENDIF 279 311 #endif … … 287 319 ! 288 320 ENDIF 289 ! 290 ! write out mesh mask 291 IF ( nn_msh > 0 ) CALL dom_wri 321 322 ! User defined alteration of fmask (use to reduce ocean transport in specified straits) 323 ! -------------------------------- 324 ! 325 326 CALL usr_def_fmask( cn_cfg, nn_cfg, fmask ) 292 327 ! 293 328 END SUBROUTINE dom_msk -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domngb.F90
r12414 r13024 11 11 !!---------------------------------------------------------------------- 12 12 USE dom_oce ! ocean space and time domain 13 USE phycst14 13 ! 15 14 USE in_out_manager ! I/O manager 16 15 USE lib_mpp ! for mppsum 16 USE phycst 17 17 18 18 IMPLICIT NONE … … 46 46 INTEGER :: ik ! working level 47 47 INTEGER , DIMENSION(2) :: iloc 48 REAL(wp) :: zlon, zmini 48 49 REAL(wp), DIMENSION(jpi,jpj) :: zglam, zgphi, zmask, zdist 49 50 !!-------------------------------------------------------------------- … … 59 60 END SELECT 60 61 61 zdist = dist(plon, plat, zglam, zgphi) 62 zlon = MOD( plon + 720., 360. ) ! plon between 0 and 360 63 zglam(:,:) = MOD( zglam(:,:) + 720., 360. ) ! glam between 0 and 360 64 IF( zlon > 270. ) zlon = zlon - 360. ! zlon between -90 and 270 65 IF( zlon < 90. ) WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360. ! glam between -180 and 180 66 zglam(:,:) = zglam(:,:) - zlon 62 67 63 68 IF( lk_mpp ) THEN -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domutl.F90
r12414 r13024 54 54 INTEGER :: iim1, ijm1, ikm1 ! working integer 55 55 INTEGER :: nseed ! size of the stack 56 TYPE (idx), POINTER :: seed 56 TYPE (idx), POINTER :: seed => NULL() 57 57 !!---------------------------------------------------------------------- 58 58 ! 59 59 ! initialisation seed 60 NULLIFY(seed)60 ! NULLIFY(seed) 61 61 ! 62 62 ! create the first seed and update nseed for all processors … … 150 150 INTEGER :: iim1, ijm1 ! working integer 151 151 INTEGER :: nseed ! size of the stack 152 TYPE (idx), POINTER :: seed 152 TYPE (idx), POINTER :: seed => NULL() 153 153 !!---------------------------------------------------------------------- 154 154 ! 155 155 ! initialisation seed 156 NULLIFY(seed)156 !NULLIFY(seed) 157 157 ! 158 158 ! create the first seed and update nseed for all processors … … 259 259 zpt_tmp => pt_idx%next 260 260 DEALLOCATE(pt_idx) 261 NULLIFY(pt_idx) 261 !NULLIFY(pt_idx) 262 pt_idx => NULL() 262 263 pt_idx => zpt_tmp 263 264 ELSE 264 NULLIFY(pt_idx) 265 !NULLIFY(pt_idx) 266 pt_idx => NULL() 265 267 ENDIF 266 268 END SUBROUTINE del_head_idx -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domzgr.F90
r12414 r13024 17 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capability 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 20 20 !! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying 21 21 !!---------------------------------------------------------------------- … … 36 36 !!--------------------------------------------------------------------- 37 37 USE dom_oce ! ocean domain 38 USE depth_e3 ! depth <=> e3 38 39 ! 39 40 USE in_out_manager ! I/O manager … … 92 93 CONTAINS 93 94 94 SUBROUTINE dom_zgr 95 SUBROUTINE dom_zgr( k_top, k_bot ) 95 96 !!---------------------------------------------------------------------- 96 97 !! *** ROUTINE dom_zgr *** … … 109 110 !! ** Action : define gdep., e3., mbathy and bathy 110 111 !!---------------------------------------------------------------------- 112 INTEGER, DIMENSION(:,:), INTENT(out) :: k_top, k_bot ! ocean first and last level indices 113 ! 111 114 INTEGER :: ioptio, ibat ! local integer 112 115 INTEGER :: ios 113 116 ! 117 INTEGER :: jk 118 REAL(wp) :: zrefdep ! depth of the reference level (~10m) 119 120 114 121 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 115 122 !!---------------------------------------------------------------------- … … 124 131 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp ) 125 132 IF(lwm) WRITE ( numond, namzgr ) 133 134 IF(ln_read_cfg) THEN 135 IF(lwp) WRITE(numout,*) 136 IF(lwp) WRITE(numout,*) ' ==>>> Read vertical mesh in ', TRIM( cn_domcfg ), ' file' 137 ! 138 CALL zgr_read ( ln_zco , ln_zps , ln_sco, ln_isfcav, & 139 & gdept_1d, gdepw_1d, e3t_1d, e3w_1d , & ! 1D gridpoints depth 140 & gdept_0 , gdepw_0 , & ! gridpoints depth 141 & e3t_0 , e3u_0 , e3v_0 , e3f_0 , & ! vertical scale factors 142 & e3w_0 , e3uw_0 , e3vw_0 , & ! vertical scale factors 143 & k_top , k_bot ) ! 1st & last ocean level 144 ! 145 !!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears 146 ! ! Compute gde3w_0 (vertical sum of e3w) 147 ! gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 148 ! DO jk = 2, jpk 149 ! gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 150 ! END DO 151 ! 152 153 ! ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top) 154 CALL zgr_top_bot( k_top, k_bot ) ! with a minimum value set to 1 155 156 ! ! deepest/shallowest W level Above/Below ~10m 157 !!gm BUG in s-coordinate this does not work! 158 zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness) 159 nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 160 nla10 = nlb10 - 1 ! deepest W level Above ~10m 161 !!gm end bug 162 163 ENDIF 126 164 127 165 IF(lwp) THEN ! Control print … … 150 188 IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1 151 189 IF( ioptio > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' ) 190 191 IF(.NOT.ln_read_cfg) THEN 152 192 ! 153 193 ! Build the vertical coordinate system … … 163 203 ! ----------------------------------- 164 204 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 205 k_bot = mbkt 165 206 CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points 207 k_top = mikt 208 209 ENDIF 166 210 ! 167 211 IF( nprint == 1 .AND. lwp ) THEN 168 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 212 WRITE(numout,*) ' MIN val k_top ', MINVAL( k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) ) 213 WRITE(numout,*) ' MIN val k_bot ', MINVAL( k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) ) 169 214 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 170 215 & ' w ', MINVAL( gdepw_0(:,:,:) ) … … 181 226 & ' w ', MAXVAL( e3w_0(:,:,:) ) 182 227 ENDIF 228 183 229 ! 184 230 END SUBROUTINE dom_zgr 185 231 232 SUBROUTINE zgr_read( ld_zco , ld_zps , ld_sco , ld_isfcav, & ! type of vertical coordinate 233 & pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d , & ! 1D reference vertical coordinate 234 & pdept , pdepw , & ! 3D t & w-points depth 235 & pe3t , pe3u , pe3v , pe3f , & ! vertical scale factors 236 & pe3w , pe3uw , pe3vw , & ! - - - 237 & k_top , k_bot ) ! top & bottom ocean level 238 !!--------------------------------------------------------------------- 239 !! *** ROUTINE zgr_read *** 240 !! 241 !! ** Purpose : Read the vertical information in the domain configuration file 242 !! 243 !!---------------------------------------------------------------------- 244 LOGICAL , INTENT(out) :: ld_zco, ld_zps, ld_sco ! vertical coordinate flags 245 LOGICAL , INTENT(out) :: ld_isfcav ! under iceshelf cavity flag 246 REAL(wp), DIMENSION(:) , INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m] 247 REAL(wp), DIMENSION(:) , INTENT(out) :: pe3t_1d , pe3w_1d ! 1D vertical scale factors [m] 248 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pdept, pdepw ! grid-point depth [m] 249 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3t , pe3u , pe3v , pe3f ! vertical scale factors [m] 250 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3w , pe3uw, pe3vw ! - - - 251 INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top , k_bot ! first & last ocean level 252 ! 253 INTEGER :: jk ! dummy loop index 254 INTEGER :: inum ! local logical unit 255 REAL(WP) :: z_zco, z_zps, z_sco, z_cav 256 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 257 !!---------------------------------------------------------------------- 258 ! 259 IF(lwp) THEN 260 WRITE(numout,*) 261 WRITE(numout,*) ' zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file' 262 WRITE(numout,*) ' ~~~~~~~~' 263 ENDIF 264 ! 265 CALL iom_open( cn_domcfg, inum ) 266 ! 267 ! !* type of vertical coordinate 268 CALL iom_get( inum, 'ln_zco' , z_zco ) 269 CALL iom_get( inum, 'ln_zps' , z_zps ) 270 CALL iom_get( inum, 'ln_sco' , z_sco ) 271 IF( z_zco == 0._wp ) THEN ; ld_zco = .false. ; ELSE ; ld_zco = .true. ; ENDIF 272 IF( z_zps == 0._wp ) THEN ; ld_zps = .false. ; ELSE ; ld_zps = .true. ; ENDIF 273 IF( z_sco == 0._wp ) THEN ; ld_sco = .false. ; ELSE ; ld_sco = .true. ; ENDIF 274 ! 275 ! !* ocean cavities under iceshelves 276 CALL iom_get( inum, 'ln_isfcav', z_cav ) 277 IF( z_cav == 0._wp ) THEN ; ld_isfcav = .false. ; ELSE ; ld_isfcav = .true. ; ENDIF 278 ! 279 ! !* vertical scale factors 280 CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , pe3t_1d ) ! 1D reference coordinate 281 CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , pe3w_1d ) 282 ! 283 CALL iom_get( inum, jpdom_data, 'e3t_0' , pe3t , lrowattr=ln_use_jattr ) ! 3D coordinate 284 CALL iom_get( inum, jpdom_data, 'e3u_0' , pe3u , lrowattr=ln_use_jattr ) 285 CALL iom_get( inum, jpdom_data, 'e3v_0' , pe3v , lrowattr=ln_use_jattr ) 286 CALL iom_get( inum, jpdom_data, 'e3f_0' , pe3f , lrowattr=ln_use_jattr ) 287 CALL iom_get( inum, jpdom_data, 'e3w_0' , pe3w , lrowattr=ln_use_jattr ) 288 CALL iom_get( inum, jpdom_data, 'e3uw_0' , pe3uw , lrowattr=ln_use_jattr ) 289 CALL iom_get( inum, jpdom_data, 'e3vw_0' , pe3vw , lrowattr=ln_use_jattr ) 290 ! 291 ! !* depths 292 ! !- old depth definition (obsolescent feature) 293 IF( iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0 .AND. & 294 & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0 .AND. & 295 & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0 .AND. & 296 & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0 ) THEN 297 CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', & 298 & ' depths at t- and w-points read in the domain configuration file') 299 CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d ) 300 CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d ) 301 CALL iom_get( inum, jpdom_data , 'gdept_0' , pdept , lrowattr=ln_use_jattr ) 302 CALL iom_get( inum, jpdom_data , 'gdepw_0' , pdepw , lrowattr=ln_use_jattr ) 303 ! 304 ELSE !- depths computed from e3. scale factors 305 CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d ) ! 1D reference depth 306 CALL e3_to_depth( pe3t , pe3w , pdept , pdepw ) ! 3D depths 307 IF(lwp) THEN 308 WRITE(numout,*) 309 WRITE(numout,*) ' Reference 1D z-coordinate depth and scale factors:' 310 WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" ) 311 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk ) 312 ENDIF 313 ENDIF 314 ! 315 ! !* ocean top and bottom level 316 CALL iom_get( inum, jpdom_data, 'top_level' , z2d , lrowattr=ln_use_jattr ) ! 1st wet T-points (ISF) 317 k_top(:,:) = NINT( z2d(:,:) ) 318 CALL iom_get( inum, jpdom_data, 'bottom_level' , z2d , lrowattr=ln_use_jattr ) ! last wet T-points 319 k_bot(:,:) = NINT( z2d(:,:) ) 320 ! 321 ! reference depth for negative bathy (wetting and drying only) 322 ! IF( ll_wd ) CALL iom_get( inum, 'rn_wd_ref_depth' , ssh_ref ) 323 ! 324 CALL iom_close( inum ) 325 ! 326 END SUBROUTINE zgr_read 327 328 SUBROUTINE zgr_top_bot( k_top, k_bot ) 329 !!---------------------------------------------------------------------- 330 !! *** ROUTINE zgr_top_bot *** 331 !! 332 !! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays) 333 !! 334 !! ** Method : computes from k_top and k_bot with a minimum value of 1 over land 335 !! 336 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 337 !! ocean level at t-, u- & v-points 338 !! (min value = 1) 339 !! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest 340 !! ocean level at t-, u- & v-points 341 !! (min value = 1 over land) 342 !!---------------------------------------------------------------------- 343 INTEGER , DIMENSION(:,:), INTENT(in) :: k_top, k_bot ! top & bottom ocean level indices 344 ! 345 INTEGER :: ji, jj ! dummy loop indices 346 REAL(wp), DIMENSION(jpi,jpj) :: zk ! workspace 347 !!---------------------------------------------------------------------- 348 ! 349 IF(lwp) WRITE(numout,*) 350 IF(lwp) WRITE(numout,*) ' zgr_top_bot : ocean top and bottom k-index of T-, U-, V- and W-levels ' 351 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~' 352 ! 353 mikt(:,:) = MAX( k_top(:,:) , 1 ) ! top ocean k-index of T-level (=1 over land) 354 ! 355 mbkt(:,:) = MAX( k_bot(:,:) , 1 ) ! bottom ocean k-index of T-level (=1 over land) 356 357 ! ! N.B. top k-index of W-level = mikt 358 ! ! bottom k-index of W-level = mbkt+1 359 DO jj = 1, jpjm1 360 DO ji = 1, jpim1 361 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 362 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 363 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 364 ! 365 mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) ) 366 mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) ) 367 END DO 368 END DO 369 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 370 zk(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'U', 1. ) ; miku(:,:) = MAX( NINT( zk(:,:) ), 1 ) 371 zk(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'V', 1. ) ; mikv(:,:) = MAX( NINT( zk(:,:) ), 1 ) 372 zk(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'F', 1. ) ; mikf(:,:) = MAX( NINT( zk(:,:) ), 1 ) 373 ! 374 zk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'U', 1. ) ; mbku(:,:) = MAX( NINT( zk(:,:) ), 1 ) 375 zk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'V', 1. ) ; mbkv(:,:) = MAX( NINT( zk(:,:) ), 1 ) 376 ! 377 END SUBROUTINE zgr_top_bot 186 378 187 379 SUBROUTINE zgr_z … … 665 857 ENDIF 666 858 667 zbathy(:,:) = FLOAT( mbathy(:,:) ) 668 CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp ) 669 mbathy(:,:) = INT( zbathy(:,:) ) 670 859 IF( lk_mpp ) THEN 860 zbathy(:,:) = FLOAT( mbathy(:,:) ) 861 CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp ) 862 mbathy(:,:) = INT( zbathy(:,:) ) 863 ENDIF 671 864 ! ! East-west cyclic boundary conditions 672 865 IF( jperio == 0 ) THEN … … 1061 1254 END SUBROUTINE zgr_zps 1062 1255 1256 1063 1257 SUBROUTINE zgr_sco 1064 1258 !!---------------------------------------------------------------------- … … 1262 1456 END DO 1263 1457 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1264 CALL lbc_lnk( ' domzgr',zenv, 'T', 1._wp, 'no0' )1458 CALL lbc_lnk( 'toto',zenv, 'T', 1._wp, 'no0' ) 1265 1459 ! ! ================ ! 1266 1460 END DO ! End loop ! -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/iom_nf90.F90
r12414 r13024 530 530 INTEGER :: idim3 ! id of the third dimension 531 531 ! 532 INTEGER :: nldi_save, nlei_save !:patch before we remove periodicity and close boundaries in output files533 INTEGER :: nldj_save, nlej_save !:532 ! INTEGER :: nldi_save, nlei_save !:patch before we remove periodicity and close boundaries in output files 533 ! INTEGER :: nldj_save, nlej_save !: 534 534 !--------------------------------------------------------------------- 535 535 ! … … 540 540 ! without this, issue in some model decomposition 541 541 ! seb: patch before we remove periodicity and close boundaries in output files 542 nldi_save = nldi ; nlei_save = nlei543 nldj_save = nldj ; nlej_save = nlej544 IF( nimpp == 1 ) nldi = 1545 IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi546 IF( njmpp == 1 ) nldj = 1547 IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj542 ! nldi_save = nldi ; nlei_save = nlei 543 ! nldj_save = nldj ; nlej_save = nlej 544 ! IF( nimpp == 1 ) nldi = 1 545 ! IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi 546 ! IF( njmpp == 1 ) nldj = 1 547 ! IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj 548 548 ! 549 549 ! define dimension variables if it is not already done … … 719 719 ENDIF 720 720 ! 721 nldi = nldi_save ; nlei = nlei_save722 nldj = nldj_save ; nlej = nlej_save721 ! nldi = nldi_save ; nlei = nlei_save 722 ! nldj = nldj_save ; nlej = nlej_save 723 723 ! 724 724 END SUBROUTINE iom_nf90_rp0123d -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/lib_mpp.F90
r12414 r13024 1662 1662 CHARACTER(len=*),DIMENSION(:) :: ldtxt 1663 1663 CHARACTER(len=*) :: ldname 1664 INTEGER :: function_value1665 1664 INTEGER :: kumnam_ref, knumnam_cfg , kumond , kstop 1666 1665 IF( PRESENT( localComm ) ) mpi_comm_oce = localComm -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/mpp_loc_generic.h90
r12414 r13024 17 17 # define MPI_OPERATION mpi_maxloc 18 18 # define LOC_OPERATION MAXLOC 19 # define ERRVAL -HUGE20 19 # endif 21 20 # if defined OPERATION_MINLOC 22 21 # define MPI_OPERATION mpi_minloc 23 22 # define LOC_OPERATION MINLOC 24 # define ERRVAL HUGE25 23 # endif 26 24 … … 44 42 ! 45 43 idim = SIZE(kindex) 44 ALLOCATE ( ilocs(idim) ) 46 45 ! 47 IF ( ALL(MASK_IN(:,:,:) /= 1._wp) ) THEN 48 ! special case for land processors 49 zmin = ERRVAL(zmin) 50 index0 = 0 51 ELSE 52 ALLOCATE ( ilocs(idim) ) 53 ! 54 ilocs = LOC_OPERATION( ARRAY_IN(:,:,:) , mask= MASK_IN(:,:,:) == 1._wp ) 55 zmin = ARRAY_IN(ilocs(1),ilocs(2),ilocs(3)) 56 ! 57 kindex(1) = mig( ilocs(1) ) 46 ilocs = LOC_OPERATION( ARRAY_IN(:,:,:) , mask= MASK_IN(:,:,:) == 1._wp ) 47 zmin = ARRAY_IN(ilocs(1),ilocs(2),ilocs(3)) 48 ! 49 kindex(1) = ilocs(1) + nimpp - 1 58 50 # if defined DIM_2d || defined DIM_3d /* avoid warning when kindex has 1 element */ 59 kindex(2) = mjg( ilocs(2) )51 kindex(2) = ilocs(2) + njmpp - 1 60 52 # endif 61 53 # if defined DIM_3d /* avoid warning when kindex has 2 elements */ 62 54 kindex(3) = ilocs(3) 63 55 # endif 64 65 66 67 56 ! 57 DEALLOCATE (ilocs) 58 ! 59 index0 = kindex(1)-1 ! 1d index starting at 0 68 60 # if defined DIM_2d || defined DIM_3d /* avoid warning when kindex has 1 element */ 69 61 index0 = index0 + jpiglo * (kindex(2)-1) 70 62 # endif 71 63 # if defined DIM_3d /* avoid warning when kindex has 2 elements */ 72 64 index0 = index0 + jpiglo * jpjglo * (kindex(3)-1) 73 65 # endif 74 END IF75 66 zain(1,:) = zmin 76 67 zain(2,:) = REAL(index0, wp) … … 107 98 #undef LOC_OPERATION 108 99 #undef INDEX_TYPE 109 #undef ERRVAL -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/mppini.F90
r12414 r13024 99 99 CALL ctl_stop( 'mpp_init: equality jpni = jpnj = jpnij = 1 is not satisfied', & 100 100 & 'the domain is lay out for distributed memory computing!' ) 101 102 #if defined key_agrif 103 IF (.not.agrif_root()) THEN 104 CALL agrif_nemo_init 105 ENDIF 106 107 IF( .NOT. Agrif_Root() ) THEN ! AGRIF children: specific setting (cf. agrif_user.F90) 108 print *,'nbcellsx = ',nbcellsx,nbghostcells_x 109 print *,'nbcellsy = ',nbcellsy,nbghostcells_y_s,nbghostcells_y_n 110 IF( jpiglo /= nbcellsx + 2 + 2*nbghostcells_x ) THEN 111 IF(lwp) THEN 112 WRITE(numout,*) 113 WRITE(numout,*) 'jpiglo should be: ', nbcellsx + 2 + 2*nbghostcells_x 114 ENDIF 115 CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells_x' ) 116 ENDIF 117 IF( jpjglo /= nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n ) THEN 118 IF(lwp) THEN 119 WRITE(numout,*) 120 WRITE(numout,*) 'jpjglo shoud be: ', nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n 121 ENDIF 122 CALL ctl_stop( 'STOP', & 123 'mpp_init: Agrif children requires jpjglo == nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n' ) 124 ENDIF 125 IF( ln_use_jattr ) CALL ctl_stop( 'STOP', 'mpp_init:Agrif children requires ln_use_jattr = .false. ' ) 126 ENDIF 127 #endif 101 128 ! 102 129 END SUBROUTINE mpp_init … … 274 301 IF(lwp) THEN 275 302 WRITE(numout,*) 276 WRITE(numout,*) 'jpiglo shou d be: ', nbcellsx + 2 + 2*nbghostcells303 WRITE(numout,*) 'jpiglo should be: ', nbcellsx + 2 + 2*nbghostcells 277 304 ENDIF 278 305 CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells' ) -
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/nemogcm.F90
r12414 r13024 408 408 !!---------------------------------------------------------------------- 409 409 ! 410 ierr = dom_oce_alloc () ! ocean domain 410 ierr = 0 411 ierr = ierr + dom_oce_alloc () ! ocean domain 411 412 ! 412 413 CALL mpp_sum( 'nemogcm', ierr ) … … 416 417 417 418 419 SUBROUTINE nemo_partition( num_pes ) 420 !!---------------------------------------------------------------------- 421 !! *** ROUTINE nemo_partition *** 422 !! 423 !! ** Purpose : 424 !! 425 !! ** Method : 426 !!---------------------------------------------------------------------- 427 INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have 428 ! 429 INTEGER, PARAMETER :: nfactmax = 20 430 INTEGER :: nfact ! The no. of factors returned 431 INTEGER :: ierr ! Error flag 432 INTEGER :: ji 433 INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are closest in value 434 INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors 435 !!---------------------------------------------------------------------- 436 ! 437 ierr = 0 438 ! 439 CALL factorise( ifact, nfactmax, nfact, num_pes, ierr ) 440 ! 441 IF( nfact <= 1 ) THEN 442 WRITE (numout, *) 'WARNING: factorisation of number of PEs failed' 443 WRITE (numout, *) ' : using grid of ',num_pes,' x 1' 444 jpnj = 1 445 jpni = num_pes 446 ELSE 447 ! Search through factors for the pair that are closest in value 448 mindiff = 1000000 449 imin = 1 450 DO ji = 1, nfact-1, 2 451 idiff = ABS( ifact(ji) - ifact(ji+1) ) 452 IF( idiff < mindiff ) THEN 453 mindiff = idiff 454 imin = ji 455 ENDIF 456 END DO 457 jpnj = ifact(imin) 458 jpni = ifact(imin + 1) 459 ENDIF 460 ! 461 jpnij = jpni*jpnj 462 ! 463 END SUBROUTINE nemo_partition 464 465 466 SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr ) 467 !!---------------------------------------------------------------------- 468 !! *** ROUTINE factorise *** 469 !! 470 !! ** Purpose : return the prime factors of n. 471 !! knfax factors are returned in array kfax which is of 472 !! maximum dimension kmaxfax. 473 !! ** Method : 474 !!---------------------------------------------------------------------- 475 INTEGER , INTENT(in ) :: kn, kmaxfax 476 INTEGER , INTENT( out) :: kerr, knfax 477 INTEGER, DIMENSION(kmaxfax), INTENT( out) :: kfax 478 ! 479 INTEGER :: ifac, jl, inu 480 INTEGER, PARAMETER :: ntest = 14 481 INTEGER, DIMENSION(ntest) :: ilfax 482 !!---------------------------------------------------------------------- 483 ! 484 ! lfax contains the set of allowed factors. 485 ilfax(:) = (/(2**jl,jl=ntest,1,-1)/) 486 ! 487 ! Clear the error flag and initialise output vars 488 kerr = 0 489 kfax = 1 490 knfax = 0 491 ! 492 ! Find the factors of n. 493 IF( kn == 1 ) GOTO 20 494 495 ! nu holds the unfactorised part of the number. 496 ! knfax holds the number of factors found. 497 ! l points to the allowed factor list. 498 ! ifac holds the current factor. 499 ! 500 inu = kn 501 knfax = 0 502 ! 503 DO jl = ntest, 1, -1 504 ! 505 ifac = ilfax(jl) 506 IF( ifac > inu ) CYCLE 507 508 ! Test whether the factor will divide. 509 510 IF( MOD(inu,ifac) == 0 ) THEN 511 ! 512 knfax = knfax + 1 ! Add the factor to the list 513 IF( knfax > kmaxfax ) THEN 514 kerr = 6 515 write (*,*) 'FACTOR: insufficient space in factor array ', knfax 516 return 517 ENDIF 518 kfax(knfax) = ifac 519 ! Store the other factor that goes with this one 520 knfax = knfax + 1 521 kfax(knfax) = inu / ifac 522 !WRITE (*,*) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax) 523 ENDIF 524 ! 525 END DO 526 ! 527 20 CONTINUE ! Label 20 is the exit point from the factor search loop. 528 ! 529 END SUBROUTINE factorise 530 531 532 SUBROUTINE nemo_northcomms 533 !!---------------------------------------------------------------------- 534 !! *** ROUTINE nemo_northcomms *** 535 !! ** Purpose : Setup for north fold exchanges with explicit 536 !! point-to-point messaging 537 !! 538 !! ** Method : Initialization of the northern neighbours lists. 539 !!---------------------------------------------------------------------- 540 !! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE) 541 !! 2.0 ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC) 542 !!---------------------------------------------------------------------- 543 INTEGER :: sxM, dxM, sxT, dxT, jn 544 INTEGER :: njmppmax 545 !!---------------------------------------------------------------------- 546 ! 547 njmppmax = MAXVAL( njmppt ) 548 ! 549 !initializes the north-fold communication variables 550 isendto(:) = 0 551 nsndto = 0 552 ! 553 !if I am a process in the north 554 IF ( njmpp == njmppmax ) THEN 555 !sxM is the first point (in the global domain) needed to compute the 556 !north-fold for the current process 557 sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1 558 !dxM is the last point (in the global domain) needed to compute the 559 !north-fold for the current process 560 dxM = jpiglo - nimppt(narea) + 2 561 562 !loop over the other north-fold processes to find the processes 563 !managing the points belonging to the sxT-dxT range 564 565 DO jn = 1, jpni 566 !sxT is the first point (in the global domain) of the jn 567 !process 568 sxT = nfiimpp(jn, jpnj) 569 !dxT is the last point (in the global domain) of the jn 570 !process 571 dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1 572 IF ((sxM .gt. sxT) .AND. (sxM .lt. dxT)) THEN 573 nsndto = nsndto + 1 574 isendto(nsndto) = jn 575 ELSEIF ((sxM .le. sxT) .AND. (dxM .ge. dxT)) THEN 576 nsndto = nsndto + 1 577 isendto(nsndto) = jn 578 ELSEIF ((dxM .lt. dxT) .AND. (sxT .lt. dxM)) THEN 579 nsndto = nsndto + 1 580 isendto(nsndto) = jn 581 END IF