- Timestamp:
- 2020-06-26T10:26:32+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement
- Files:
-
- 46 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette@ HEADsette10 ^/utils/CI/sette@12931 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/ASM/asminc.F90
r12489 r13159 896 896 IF ( kt == nitdin_r ) THEN 897 897 ! 898 l_1st_euler = 0! Force Euler forward step898 l_1st_euler = .TRUE. ! Force Euler forward step 899 899 ! 900 900 ! Sea-ice : SI3 case -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/BDY/bdydta.F90
r12547 r13159 91 91 INTEGER :: jbdy, jfld, jstart, jend, ib, jl ! dummy loop indices 92 92 INTEGER :: ii, ij, ik, igrd, ipl ! local integers 93 INTEGER, DIMENSION(jpbgrd) :: ilen194 93 TYPE(OBC_DATA) , POINTER :: dta_alias ! short cut 95 94 TYPE(FLD), DIMENSION(:), POINTER :: bf_alias … … 116 115 END DO 117 116 ENDIF 118 IF( dta_bdy(jbdy)%lneed_dyn2d) THEN117 IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain 119 118 igrd = 2 120 DO ib = 1, SIZE(dta_bdy(jbdy)%u2d) ! u2d is used only on the rim except if ln_full_vel = T, see bdy_dta_init119 DO ib = 1, SIZE(dta_bdy(jbdy)%u2d) ! u2d is used either over the whole bdy or only on the rim 121 120 ii = idx_bdy(jbdy)%nbi(ib,igrd) 122 121 ij = idx_bdy(jbdy)%nbj(ib,igrd) 123 122 dta_bdy(jbdy)%u2d(ib) = uu_b(ii,ij,Kmm) * umask(ii,ij,1) 124 123 END DO 124 ENDIF 125 IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) THEN ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain 125 126 igrd = 3 126 DO ib = 1, SIZE(dta_bdy(jbdy)%v2d) ! v2d is used only on the rim except if ln_full_vel = T, see bdy_dta_init127 DO ib = 1, SIZE(dta_bdy(jbdy)%v2d) ! v2d is used either over the whole bdy or only on the rim 127 128 ii = idx_bdy(jbdy)%nbi(ib,igrd) 128 129 ij = idx_bdy(jbdy)%nbj(ib,igrd) … … 210 211 ! 211 212 ! if runoff condition: change river flow we read (in m3/s) into barotropic velocity (m/s) 212 IF( cn_tra(jbdy) == 'runoff' .AND. TRIM(bf_alias(jp_bdyu2d)%clrootname) /= 'NOT USED' ) THEN ! runoff and we read u/v2d213 IF( cn_tra(jbdy) == 'runoff' ) THEN ! runoff 213 214 ! 214 igrd = 2 ! zonal flow (m3/s) to barotropic zonal velocity (m/s) 215 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 216 ii = idx_bdy(jbdy)%nbi(ib,igrd) 217 ij = idx_bdy(jbdy)%nbj(ib,igrd) 218 dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 219 END DO 220 igrd = 3 ! meridional flow (m3/s) to barotropic meridional velocity (m/s) 221 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 222 ii = idx_bdy(jbdy)%nbi(ib,igrd) 223 ij = idx_bdy(jbdy)%nbj(ib,igrd) 224 dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 225 END DO 215 IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain 216 igrd = 2 ! zonal flow (m3/s) to barotropic zonal velocity (m/s) 217 DO ib = 1, SIZE(dta_alias%u2d) ! u2d is used either over the whole bdy or only on the rim 218 ii = idx_bdy(jbdy)%nbi(ib,igrd) 219 ij = idx_bdy(jbdy)%nbj(ib,igrd) 220 dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 221 END DO 222 ENDIF 223 IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) THEN ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain 224 igrd = 3 ! meridional flow (m3/s) to barotropic meridional velocity (m/s) 225 DO ib = 1, SIZE(dta_alias%v2d) ! v2d is used either over the whole bdy or only on the rim 226 ii = idx_bdy(jbdy)%nbi(ib,igrd) 227 ij = idx_bdy(jbdy)%nbj(ib,igrd) 228 dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 229 END DO 230 ENDIF 226 231 ENDIF 227 232 228 233 ! tidal harmonic forcing ONLY: initialise arrays 229 234 IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! we did not read ssh, u/v2d 230 IF( dta_alias%lneed_ssh) dta_alias%ssh(:) = 0._wp231 IF( dta_alias%lneed_dyn2d) dta_alias%u2d(:) = 0._wp232 IF( dta_alias%lneed_dyn2d) dta_alias%v2d(:) = 0._wp235 IF( ASSOCIATED(dta_alias%ssh) ) dta_alias%ssh(:) = 0._wp 236 IF( ASSOCIATED(dta_alias%u2d) ) dta_alias%u2d(:) = 0._wp 237 IF( ASSOCIATED(dta_alias%v2d) ) dta_alias%v2d(:) = 0._wp 233 238 ENDIF 234 239 … … 237 242 ! 238 243 igrd = 2 ! zonal velocity 239 dta_alias%u2d(:) = 0._wp ! compute barotrope zonal velocity and put it in u2d240 244 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 241 245 ii = idx_bdy(jbdy)%nbi(ib,igrd) 242 246 ij = idx_bdy(jbdy)%nbj(ib,igrd) 247 dta_alias%u2d(ib) = 0._wp ! compute barotrope zonal velocity and put it in u2d 243 248 DO ik = 1, jpkm1 244 249 dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik) … … 250 255 END DO 251 256 igrd = 3 ! meridional velocity 252 dta_alias%v2d(:) = 0._wp ! compute barotrope meridional velocity and put it in v2d253 257 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 254 258 ii = idx_bdy(jbdy)%nbi(ib,igrd) 255 259 ij = idx_bdy(jbdy)%nbj(ib,igrd) 260 dta_alias%v2d(ib) = 0._wp ! compute barotrope meridional velocity and put it in v2d 256 261 DO ik = 1, jpkm1 257 262 dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik) … … 275 280 276 281 #if defined key_si3 277 IF( dta_alias%lneed_ice ) THEN282 IF( dta_alias%lneed_ice .AND. idx_bdy(jbdy)%nblen(1) > 0 ) THEN 278 283 ! fill temperature and salinity arrays 279 284 IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) bf_alias(jp_bdyt_i)%fnow(:,1,:) = rice_tem (jbdy) … … 330 335 DO jbdy = 1, nb_bdy ! Tidal component added in ts loop 331 336 IF ( nn_dyn2d_dta(jbdy) .GE. 2 ) THEN 332 IF( cn_dyn2d(jbdy) == 'frs' ) THEN ; ilen1(:)=idx_bdy(jbdy)%nblen(:) 333 ELSE ; ilen1(:)=idx_bdy(jbdy)%nblenrim(:) 334 ENDIF 335 IF ( dta_bdy(jbdy)%lneed_ssh ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1)) 336 IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2)) 337 IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3)) 337 IF( ASSOCIATED(dta_bdy(jbdy)%ssh) ) dta_bdy_s(jbdy)%ssh(:) = dta_bdy(jbdy)%ssh(:) 338 IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) dta_bdy_s(jbdy)%u2d(:) = dta_bdy(jbdy)%u2d(:) 339 IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) dta_bdy_s(jbdy)%v2d(:) = dta_bdy(jbdy)%v2d(:) 338 340 ENDIF 339 341 END DO 340 342 ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 341 343 ! 342 ! BDY: use pt_offset=1.0 as applied at the end of the step and bdy_dta_tides is referenced at the middle of the step343 344 CALL bdy_dta_tides( kt=kt, pt_offset = 1._wp ) 344 345 ENDIF … … 348 349 ! 349 350 END SUBROUTINE bdy_dta 350 351 351 352 352 353 SUBROUTINE bdy_dta_init … … 380 381 LOGICAL :: llneed ! 381 382 LOGICAL :: llread ! 383 LOGICAL :: llfullbdy ! 382 384 TYPE(FLD_N), DIMENSION(1), TARGET :: bn_tem, bn_sal, bn_u3d, bn_v3d ! must be an array to be used with fld_fill 383 385 TYPE(FLD_N), DIMENSION(1), TARGET :: bn_ssh, bn_u2d, bn_v2d ! informations about the fields to be read … … 494 496 igrd = 2 ! U point 495 497 ipk = 1 ! surface data 496 llneed = dta_bdy(jbdy)%lneed_dyn2d ! dta_bdy(jbdy)% sshwill be needed498 llneed = dta_bdy(jbdy)%lneed_dyn2d ! dta_bdy(jbdy)%u2d will be needed 497 499 llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1 ! don't get u2d from u3d and read NetCDF file 498 500 bf_alias => bf(jp_bdyu2d,jbdy:jbdy) ! alias for u2d structure of bdy number jbdy 499 501 bn_alias => bn_u2d ! alias for u2d structure of nambdy_dta 500 IF( ln_full_vel ) THEN ; iszdim = idx_bdy(jbdy)%nblen(igrd) ! will be computed from u3d -> need on the full bdy 501 ELSE ; iszdim = idx_bdy(jbdy)%nblenrim(igrd) ! used only on the rim 502 llfullbdy = ln_full_vel .OR. cn_dyn2d(jbdy) == 'frs' ! need u2d over the whole bdy or only over the rim? 503 IF( llfullbdy ) THEN ; iszdim = idx_bdy(jbdy)%nblen(igrd) 504 ELSE ; iszdim = idx_bdy(jbdy)%nblenrim(igrd) 502 505 ENDIF 503 506 ENDIF … … 506 509 igrd = 3 ! V point 507 510 ipk = 1 ! surface data 508 llneed = dta_bdy(jbdy)%lneed_dyn2d ! dta_bdy(jbdy)% sshwill be needed511 llneed = dta_bdy(jbdy)%lneed_dyn2d ! dta_bdy(jbdy)%v2d will be needed 509 512 llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1 ! don't get v2d from v3d and read NetCDF file 510 513 bf_alias => bf(jp_bdyv2d,jbdy:jbdy) ! alias for v2d structure of bdy number jbdy 511 514 bn_alias => bn_v2d ! alias for v2d structure of nambdy_dta 512 IF( ln_full_vel ) THEN ; iszdim = idx_bdy(jbdy)%nblen(igrd) ! will be computed from v3d -> need on the full bdy 513 ELSE ; iszdim = idx_bdy(jbdy)%nblenrim(igrd) ! used only on the rim 515 llfullbdy = ln_full_vel .OR. cn_dyn2d(jbdy) == 'frs' ! need v2d over the whole bdy or only over the rim? 516 IF( llfullbdy ) THEN ; iszdim = idx_bdy(jbdy)%nblen(igrd) 517 ELSE ; iszdim = idx_bdy(jbdy)%nblenrim(igrd) 514 518 ENDIF 515 519 ENDIF -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/BDY/bdyini.F90
r12377 r13159 19 19 USE oce ! ocean dynamics and tracers variables 20 20 USE dom_oce ! ocean space and time domain 21 USE sbc_oce , ONLY: nn_ice 21 22 USE bdy_oce ! unstructured open boundary conditions 22 23 USE bdydta ! open boundary cond. setting (bdy_dta_init routine) 23 24 USE bdytides ! open boundary cond. setting (bdytide_init routine) 24 25 USE tide_mod, ONLY: ln_tide ! tidal forcing 25 USE phycst 26 USE phycst , ONLY: rday 26 27 ! 27 28 USE in_out_manager ! I/O units … … 315 316 316 317 dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none' 318 319 IF( dta_bdy(ib_bdy)%lneed_ice .AND. nn_ice /= 2 ) THEN 320 WRITE(ctmp1,*) 'bdy number ', ib_bdy,', needs ice model but nn_ice = ', nn_ice 321 CALL ctl_stop( ctmp1 ) 322 ENDIF 317 323 318 324 IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/BDY/bdytides.F90
r12489 r13159 65 65 !! namelist variables 66 66 !!------------------- 67 CHARACTER(len=80) :: filtide ! :Filename root for tidal input files68 LOGICAL :: ln_bdytide_2ddta ! :If true, read 2d harmonic data67 CHARACTER(len=80) :: filtide ! Filename root for tidal input files 68 LOGICAL :: ln_bdytide_2ddta ! If true, read 2d harmonic data 69 69 !! 70 INTEGER :: ib_bdy, itide, ib ! :dummy loop indices71 INTEGER :: ii, ij ! :dummy loop indices70 INTEGER :: ib_bdy, itide, ib ! dummy loop indices 71 INTEGER :: ii, ij ! dummy loop indices 72 72 INTEGER :: inum, igrd 73 INTEGER , DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays)73 INTEGER :: isz ! bdy data size 74 74 INTEGER :: ios ! Local integer output status for namelist read 75 75 INTEGER :: nbdy_rdstart, nbdy_loc 76 CHARACTER(LEN=50) :: cerrmsg ! :error string77 CHARACTER(len=80) :: clfile ! :full file name for tidal input file78 REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: dta_read ! :work space to read in tidal harmonics data79 REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: ztr, zti ! :" " " " " " " "76 CHARACTER(LEN=50) :: cerrmsg ! error string 77 CHARACTER(len=80) :: clfile ! full file name for tidal input file 78 REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: dta_read ! work space to read in tidal harmonics data 79 REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: ztr, zti ! " " " " " " " " 80 80 !! 81 TYPE(TIDES_DATA), POINTER :: td !: local short cut 81 TYPE(TIDES_DATA), POINTER :: td ! local short cut 82 TYPE( OBC_DATA), POINTER :: dta ! local short cut 82 83 !! 83 84 NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta … … 93 94 IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN 94 95 ! 95 td => tides(ib_bdy) 96 96 td => tides(ib_bdy) 97 dta => dta_bdy(ib_bdy) 98 97 99 ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 98 100 filtide(:) = '' … … 130 132 IF(lwp) WRITE(numout,*) ' ' 131 133 132 ! Allocate space for tidal harmonics data - get size from OBC data arrays 134 ! Allocate space for tidal harmonics data - get size from BDY data arrays 135 ! Allocate also slow varying data in the case of time splitting: 136 ! Do it anyway because at this stage knowledge of free surface scheme is unknown 133 137 ! ----------------------------------------------------------------------- 134 135 ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 136 ! relaxation area 137 IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN ; ilen0(:) = idx_bdy(ib_bdy)%nblen (:) 138 ELSE ; ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:) 139 ENDIF 140 141 ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) ) 142 ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) ) 143 144 ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) ) 145 ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) ) 146 147 ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) ) 148 ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 149 150 td%ssh0(:,:,:) = 0._wp 151 td%ssh (:,:,:) = 0._wp 152 td%u0 (:,:,:) = 0._wp 153 td%u (:,:,:) = 0._wp 154 td%v0 (:,:,:) = 0._wp 155 td%v (:,:,:) = 0._wp 156 138 IF( ASSOCIATED(dta%ssh) ) THEN ! we use bdy ssh on this mpi subdomain 139 isz = SIZE(dta%ssh) 140 ALLOCATE( td%ssh0( isz, nb_harmo, 2 ), td%ssh( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%ssh( isz ) ) 141 dta_bdy_s(ib_bdy)%ssh(:) = 0._wp ! needed? 142 ENDIF 143 IF( ASSOCIATED(dta%u2d) ) THEN ! we use bdy u2d on this mpi subdomain 144 isz = SIZE(dta%u2d) 145 ALLOCATE( td%u0 ( isz, nb_harmo, 2 ), td%u ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%u2d( isz ) ) 146 dta_bdy_s(ib_bdy)%u2d(:) = 0._wp ! needed? 147 ENDIF 148 IF( ASSOCIATED(dta%v2d) ) THEN ! we use bdy v2d on this mpi subdomain 149 isz = SIZE(dta%v2d) 150 ALLOCATE( td%v0 ( isz, nb_harmo, 2 ), td%v ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%v2d( isz ) ) 151 dta_bdy_s(ib_bdy)%v2d(:) = 0._wp ! needed? 152 ENDIF 153 154 ! fill td%ssh0, td%u0, td%v0 155 ! ----------------------------------------------------------------------- 157 156 IF( ln_bdytide_2ddta ) THEN 157 ! 158 158 ! It is assumed that each data file contains all complex harmonic amplitudes 159 159 ! given on the global domain (ie global, jpiglo x jpjglo) … … 162 162 ! 163 163 ! SSH fields 164 clfile = TRIM(filtide)//'_grid_T.nc' 165 CALL iom_open( clfile , inum ) 166 igrd = 1 ! Everything is at T-points here 167 DO itide = 1, nb_harmo 168 CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) ) 169 CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 170 DO ib = 1, ilen0(igrd) 171 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 172 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 173 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove? 174 td%ssh0(ib,itide,1) = ztr(ii,ij) 175 td%ssh0(ib,itide,2) = zti(ii,ij) 176 END DO 177 END DO 178 CALL iom_close( inum ) 164 IF( ASSOCIATED(dta%ssh) ) THEN ! we use bdy ssh on this mpi subdomain 165 clfile = TRIM(filtide)//'_grid_T.nc' 166 CALL iom_open( clfile , inum ) 167 igrd = 1 ! Everything is at T-points here 168 DO itide = 1, nb_harmo 169 CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) ) 170 CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 171 DO ib = 1, SIZE(dta%ssh) 172 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 173 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 174 td%ssh0(ib,itide,1) = ztr(ii,ij) 175 td%ssh0(ib,itide,2) = zti(ii,ij) 176 END DO 177 END DO 178 CALL iom_close( inum ) 179 ENDIF 179 180 ! 180 181 ! U fields 181 clfile = TRIM(filtide)//'_grid_U.nc' 182 CALL iom_open( clfile , inum ) 183 igrd = 2 ! Everything is at U-points here 184 DO itide = 1, nb_harmo 185 CALL iom_get ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) ) 186 CALL iom_get ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) ) 187 DO ib = 1, ilen0(igrd) 188 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 189 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 190 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove? 191 td%u0(ib,itide,1) = ztr(ii,ij) 192 td%u0(ib,itide,2) = zti(ii,ij) 193 END DO 194 END DO 195 CALL iom_close( inum ) 182 IF( ASSOCIATED(dta%u2d) ) THEN ! we use bdy u2d on this mpi subdomain 183 clfile = TRIM(filtide)//'_grid_U.nc' 184 CALL iom_open( clfile , inum ) 185 igrd = 2 ! Everything is at U-points here 186 DO itide = 1, nb_harmo 187 CALL iom_get ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) ) 188 CALL iom_get ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) ) 189 DO ib = 1, SIZE(dta%u2d) 190 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 191 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 192 td%u0(ib,itide,1) = ztr(ii,ij) 193 td%u0(ib,itide,2) = zti(ii,ij) 194 END DO 195 END DO 196 CALL iom_close( inum ) 197 ENDIF 196 198 ! 197 199 ! V fields 198 clfile = TRIM(filtide)//'_grid_V.nc' 199 CALL iom_open( clfile , inum ) 200 igrd = 3 ! Everything is at V-points here 201 DO itide = 1, nb_harmo 202 CALL iom_get ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) ) 203 CALL iom_get ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) ) 204 DO ib = 1, ilen0(igrd) 205 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 206 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 207 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove? 208 td%v0(ib,itide,1) = ztr(ii,ij) 209 td%v0(ib,itide,2) = zti(ii,ij) 210 END DO 211 END DO 212 CALL iom_close( inum ) 200 IF( ASSOCIATED(dta%v2d) ) THEN ! we use bdy v2d on this mpi subdomain 201 clfile = TRIM(filtide)//'_grid_V.nc' 202 CALL iom_open( clfile , inum ) 203 igrd = 3 ! Everything is at V-points here 204 DO itide = 1, nb_harmo 205 CALL iom_get ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) ) 206 CALL iom_get ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) ) 207 DO ib = 1, SIZE(dta%v2d) 208 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 209 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 210 td%v0(ib,itide,1) = ztr(ii,ij) 211 td%v0(ib,itide,2) = zti(ii,ij) 212 END DO 213 END DO 214 CALL iom_close( inum ) 215 ENDIF 213 216 ! 214 217 DEALLOCATE( ztr, zti ) … … 218 221 ! Read tidal data only on bdy segments 219 222 ! 220 ALLOCATE( dta_read( MAXVAL( ilen0(1:3)), 1, 1 ) )223 ALLOCATE( dta_read( MAXVAL( idx_bdy(ib_bdy)%nblen(:) ), 1, 1 ) ) 221 224 ! 222 225 ! Open files and read in tidal forcing data … … 225 228 DO itide = 1, nb_harmo 226 229 ! ! SSH fields 227 clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc' 228 CALL iom_open( clfile, inum ) 229 CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 230 td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 231 CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 232 td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 233 CALL iom_close( inum ) 230 IF( ASSOCIATED(dta%ssh) ) THEN ! we use bdy ssh on this mpi subdomain 231 isz = SIZE(dta%ssh) 232 clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc' 233 CALL iom_open( clfile, inum ) 234 CALL fld_map( inum, 'z1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 235 td%ssh0(:,itide,1) = dta_read(1:isz,1,1) 236 CALL fld_map( inum, 'z2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 237 td%ssh0(:,itide,2) = dta_read(1:isz,1,1) 238 CALL iom_close( inum ) 239 ENDIF 234 240 ! ! U fields 235 clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc' 236 CALL iom_open( clfile, inum ) 237 CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 238 td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 239 CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 240 td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 241 CALL iom_close( inum ) 241 IF( ASSOCIATED(dta%u2d) ) THEN ! we use bdy u2d on this mpi subdomain 242 isz = SIZE(dta%u2d) 243 clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc' 244 CALL iom_open( clfile, inum ) 245 CALL fld_map( inum, 'u1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 246 td%u0(:,itide,1) = dta_read(1:isz,1,1) 247 CALL fld_map( inum, 'u2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 248 td%u0(:,itide,2) = dta_read(1:isz,1,1) 249 CALL iom_close( inum ) 250 ENDIF 242 251 ! ! V fields 243 clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc' 244 CALL iom_open( clfile, inum ) 245 CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 246 td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 247 CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 248 td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 249 CALL iom_close( inum ) 252 IF( ASSOCIATED(dta%v2d) ) THEN ! we use bdy v2d on this mpi subdomain 253 isz = SIZE(dta%v2d) 254 clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc' 255 CALL iom_open( clfile, inum ) 256 CALL fld_map( inum, 'v1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 257 td%v0(:,itide,1) = dta_read(1:isz,1,1) 258 CALL fld_map( inum, 'v2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 259 td%v0(:,itide,2) = dta_read(1:isz,1,1) 260 CALL iom_close( inum ) 261 ENDIF 250 262 ! 251 263 END DO ! end loop on tidal components … … 254 266 ! 255 267 ENDIF ! ln_bdytide_2ddta=.true. 256 !257 ! Allocate slow varying data in the case of time splitting:258 ! Do it anyway because at this stage knowledge of free surface scheme is unknown259 ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )260 ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )261 ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )262 dta_bdy_s(ib_bdy)%ssh(:) = 0._wp263 dta_bdy_s(ib_bdy)%u2d(:) = 0._wp264 dta_bdy_s(ib_bdy)%v2d(:) = 0._wp265 268 ! 266 269 ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2 … … 283 286 ! 284 287 LOGICAL :: lk_first_btstp ! =.TRUE. if time splitting and first barotropic step 285 INTEGER :: itide, ib_bdy, ib, igrd ! loop indices 286 INTEGER, DIMENSION(jpbgrd) :: ilen0 287 INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim ! short cuts 288 INTEGER :: itide, ib_bdy, ib ! loop indices 288 289 REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist, zt_offset 289 290 !!---------------------------------------------------------------------- … … 310 311 IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN 311 312 ! 312 nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)313 nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)314 !315 IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN ; ilen0(:) = nblen (:)316 ELSE ; ilen0(:) = nblenrim(:)317 ENDIF318 !319 313 ! We refresh nodal factors every day below 320 314 ! This should be done somewhere else … … 337 331 ! If time splitting, initialize arrays from slow varying open boundary data: 338 332 IF ( PRESENT(kit) ) THEN 339 IF ( dta_bdy(ib_bdy)%lneed_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))340 IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))341 IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))333 IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) dta_bdy(ib_bdy)%ssh(:) = dta_bdy_s(ib_bdy)%ssh(:) 334 IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) dta_bdy(ib_bdy)%u2d(:) = dta_bdy_s(ib_bdy)%u2d(:) 335 IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) dta_bdy(ib_bdy)%v2d(:) = dta_bdy_s(ib_bdy)%v2d(:) 342 336 ENDIF 343 337 ! … … 349 343 z_sist = zramp * SIN( z_sarg ) 350 344 ! 351 IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN 352 igrd=1 ! SSH on tracer grid 353 DO ib = 1, ilen0(igrd) 345 IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) THEN ! SSH on tracer grid 346 DO ib = 1, SIZE(dta_bdy(ib_bdy)%ssh) 354 347 dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + & 355 348 & ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + & … … 358 351 ENDIF 359 352 ! 360 IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN 361 igrd=2 ! U grid 362 DO ib = 1, ilen0(igrd) 353 IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) THEN ! U grid 354 DO ib = 1, SIZE(dta_bdy(ib_bdy)%u2d) 363 355 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + & 364 356 & ( tides(ib_bdy)%u(ib,itide,1)*z_cost + & 365 357 & tides(ib_bdy)%u(ib,itide,2)*z_sist ) 366 358 END DO 367 igrd=3 ! V grid 368 DO ib = 1, ilen0(igrd) 359 ENDIF 360 ! 361 IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) THEN ! V grid 362 DO ib = 1, SIZE(dta_bdy(ib_bdy)%v2d) 369 363 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + & 370 364 & ( tides(ib_bdy)%v(ib,itide,1)*z_cost + & … … 372 366 END DO 373 367 ENDIF 368 ! 374 369 END DO 375 END 370 ENDIF 376 371 END DO 377 372 ! … … 386 381 TYPE(TIDES_DATA), INTENT(inout) :: td ! tidal harmonics data 387 382 ! 388 INTEGER :: itide, igrd, ib ! dummy loop indices 389 INTEGER, DIMENSION(1) :: ilen0 ! length of boundary data (from OBC arrays) 383 INTEGER :: itide, isz, ib ! dummy loop indices 390 384 REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide 391 385 !!---------------------------------------------------------------------- 392 386 ! 393 igrd=1 394 ! SSH on tracer grid. 395 ilen0(1) = SIZE(td%ssh0(:,1,1)) 396 ! 397 ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) ) 398 ! 399 DO itide = 1, nb_harmo 400 DO ib = 1, ilen0(igrd) 401 mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.) 402 phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1)) 387 IF( ASSOCIATED(td%ssh0) ) THEN ! SSH on tracer grid. 388 ! 389 isz = SIZE( td%ssh0, dim = 1 ) 390 ALLOCATE( mod_tide(isz), phi_tide(isz) ) 391 ! 392 DO itide = 1, nb_harmo 393 DO ib = 1, isz 394 mod_tide(ib)=SQRT( td%ssh0(ib,itide,1)*td%ssh0(ib,itide,1) + td%ssh0(ib,itide,2)*td%ssh0(ib,itide,2) ) 395 phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1)) 396 END DO 397 DO ib = 1, isz 398 mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 399 phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u 400 END DO 401 DO ib = 1, isz 402 td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 403 td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 404 END DO 403 405 END DO 404 DO ib = 1 , ilen0(igrd) 405 mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 406 phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u 407 ENDDO 408 DO ib = 1 , ilen0(igrd) 409 td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 410 td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 411 ENDDO 412 END DO 413 ! 414 DEALLOCATE( mod_tide, phi_tide ) 406 ! 407 DEALLOCATE( mod_tide, phi_tide ) 408 ! 409 ENDIF 415 410 ! 416 411 END SUBROUTINE tide_init_elevation … … 424 419 TYPE(TIDES_DATA), INTENT(inout) :: td ! tidal harmonics data 425 420 ! 426 INTEGER :: itide, igrd, ib ! dummy loop indices 427 INTEGER, DIMENSION(3) :: ilen0 ! length of boundary data (from OBC arrays) 421 INTEGER :: itide, isz, ib ! dummy loop indices 428 422 REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide 429 423 !!---------------------------------------------------------------------- 430 424 ! 431 ilen0(2) = SIZE(td%u0(:,1,1)) 432 ilen0(3) = SIZE(td%v0(:,1,1)) 433 ! 434 igrd=2 ! U grid. 435 ! 436 ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) ) 437 ! 438 DO itide = 1, nb_harmo 439 DO ib = 1, ilen0(igrd) 440 mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.) 441 phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1)) 425 IF( ASSOCIATED(td%u0) ) THEN ! U grid. we use bdy u2d on this mpi subdomain 426 ! 427 isz = SIZE( td%u0, dim = 1 ) 428 ALLOCATE( mod_tide(isz), phi_tide(isz) ) 429 ! 430 DO itide = 1, nb_harmo 431 DO ib = 1, isz 432 mod_tide(ib)=SQRT( td%u0(ib,itide,1)*td%u0(ib,itide,1) + td%u0(ib,itide,2)*td%u0(ib,itide,2) ) 433 phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1)) 434 END DO 435 DO ib = 1, isz 436 mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 437 phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 438 END DO 439 DO ib = 1, isz 440 td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 441 td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 442 END DO 442 443 END DO 443 DO ib = 1, ilen0(igrd) 444 mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 445 phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 446 ENDDO 447 DO ib = 1, ilen0(igrd) 448 td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 449 td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 450 ENDDO 451 END DO 452 ! 453 DEALLOCATE( mod_tide , phi_tide ) 454 ! 455 igrd=3 ! V grid. 456 ! 457 ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) ) 458 459 DO itide = 1, nb_harmo 460 DO ib = 1, ilen0(igrd) 461 mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.) 462 phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1)) 444 ! 445 DEALLOCATE( mod_tide, phi_tide ) 446 ! 447 ENDIF 448 ! 449 IF( ASSOCIATED(td%v0) ) THEN ! V grid. we use bdy u2d on this mpi subdomain 450 ! 451 isz = SIZE( td%v0, dim = 1 ) 452 ALLOCATE( mod_tide(isz), phi_tide(isz) ) 453 ! 454 DO itide = 1, nb_harmo 455 DO ib = 1, isz 456 mod_tide(ib)=SQRT( td%v0(ib,itide,1)*td%v0(ib,itide,1) + td%v0(ib,itide,2)*td%v0(ib,itide,2) ) 457 phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1)) 458 END DO 459 DO ib = 1, isz 460 mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 461 phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 462 END DO 463 DO ib = 1, isz 464 td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 465 td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 466 END DO 463 467 END DO 464 DO ib = 1, ilen0(igrd) 465 mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 466 phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 467 ENDDO 468 DO ib = 1, ilen0(igrd) 469 td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 470 td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 471 ENDDO 472 END DO 473 ! 474 DEALLOCATE( mod_tide, phi_tide ) 475 ! 476 END SUBROUTINE tide_init_velocities 468 ! 469 DEALLOCATE( mod_tide, phi_tide ) 470 ! 471 ENDIF 472 ! 473 END SUBROUTINE tide_init_velocities 477 474 478 475 !!====================================================================== -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/C1D/step_c1d.F90
r12377 r13159 27 27 PRIVATE 28 28 29 PUBLIC stp_c1d ! called by opa.F9029 PUBLIC stp_c1d ! called by nemogcm.F90 30 30 31 31 !!---------------------------------------------------------------------- … … 56 56 ! 57 57 INTEGER :: jk ! dummy loop indice 58 INTEGER :: indic ! error indicator if < 059 58 !! --------------------------------------------------------------------- 60 61 indic = 0 ! reset to no error condition62 59 IF( kstp == nit000 ) CALL iom_init( "nemo") ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 63 60 IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init) … … 88 85 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 89 86 CALL dia_wri( kstp, Nnn ) ! ocean model: outputs 90 IF( lk_diahth )CALL dia_hth( kstp, Nnn ) ! Thermocline depth (20°C)87 CALL dia_hth( kstp, Nnn ) ! Thermocline depth (20°C) 91 88 92 89 … … 111 108 CALL eos( ts(:,:,:,:,Nnn), rhd, rhop, gdept_0(:,:,:) ) ! now potential density for zdfmxl 112 109 IF( ln_zdfnpc ) CALL tra_npc( kstp, Nnn, Nrhs, ts, Naa ) ! applied non penetrative convective adjustment on (t,s) 113 CALL tra_atf( kstp, Nbb, Nnn, Nrhs, Naa, ts ) ! time filtering of "now" tracer fields 114 115 110 CALL tra_atf( kstp, Nbb, Nnn, Naa, ts ) ! time filtering of "now" tracer arrays 116 111 117 112 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 139 134 ! Control and restarts 140 135 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 141 CALL stp_ctl( kstp, Nnn , indic)136 CALL stp_ctl( kstp, Nnn ) 142 137 IF( kstp == nit000 ) CALL iom_close( numror ) ! close input ocean restart file 143 138 IF( lrst_oce ) CALL rst_write( kstp, Nbb, Nnn ) ! write output ocean restart file 144 139 ! 145 140 #if defined key_iomput 146 IF( kstp == nitend .OR. indic <0 ) CALL xios_context_finalize() ! needed for XIOS141 IF( kstp == nitend .OR. nstop > 0 ) CALL xios_context_finalize() ! needed for XIOS 147 142 ! 148 143 #endif -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DIA/diaar5.F90
r12489 r13159 32 32 REAL(wp) :: vol0 ! ocean volume (interior domain) 33 33 REAL(wp) :: area_tot ! total ocean surface (interior domain) 34 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: area ! cell surface (interior domain)35 34 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: thick0 ! ocean thickness (interior domain) 36 35 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sn0 ! initial salinity … … 54 53 !!---------------------------------------------------------------------- 55 54 ! 56 ALLOCATE( area(jpi,jpj),thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc )55 ALLOCATE( thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc ) 57 56 ! 58 57 CALL mpp_sum ( 'diaar5', dia_ar5_alloc ) … … 78 77 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zarea_ssh , zbotpres ! 2D workspace 79 78 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpe, z2d ! 2D workspace 80 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zrhd , z rhop, ztpot! 3D workspace79 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zrhd , ztpot ! 3D workspace 81 80 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn ! 4D workspace 82 81 … … 88 87 IF( l_ar5 ) THEN 89 88 ALLOCATE( zarea_ssh(jpi,jpj), zbotpres(jpi,jpj), z2d(jpi,jpj) ) 90 ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk))89 ALLOCATE( zrhd(jpi,jpj,jpk) ) 91 90 ALLOCATE( ztsn(jpi,jpj,jpk,jpts) ) 92 zarea_ssh(:,:) = area(:,:) * ssh(:,:,Kmm)93 ENDIF 94 ! 95 CALL iom_put( 'e2u' , e2u (:,:) )96 CALL iom_put( 'e1v' , e1v (:,:) )97 CALL iom_put( 'areacello', area(:,:) )91 zarea_ssh(:,:) = e1e2t(:,:) * ssh(:,:,Kmm) 92 ENDIF 93 ! 94 CALL iom_put( 'e2u' , e2u (:,:) ) 95 CALL iom_put( 'e1v' , e1v (:,:) ) 96 CALL iom_put( 'areacello', e1e2t(:,:) ) 98 97 ! 99 98 IF( iom_use( 'volcello' ) .OR. iom_use( 'masscello' ) ) THEN 100 99 zrhd(:,:,jpk) = 0._wp ! ocean volume ; rhd is used as workspace 101 100 DO jk = 1, jpkm1 102 zrhd(:,:,jk) = area(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)101 zrhd(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) 103 102 END DO 104 103 CALL iom_put( 'volcello' , zrhd(:,:,:) ) ! WARNING not consistent with CMIP DR where volcello is at ca. 2000 … … 151 150 END IF 152 151 ! 153 zarho = glob_sum( 'diaar5', area(:,:) * zbotpres(:,:) )152 zarho = glob_sum( 'diaar5', e1e2t(:,:) * zbotpres(:,:) ) 154 153 zssh_steric = - zarho / area_tot 155 154 CALL iom_put( 'sshthster', zssh_steric ) 156 155 157 156 ! ! steric sea surface height 158 CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, gdept(:,:,:,Kmm) ) ! now in situ and potential density159 zrhop(:,:,jpk) = 0._wp160 CALL iom_put( 'rhop', zrhop )161 !162 157 zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice 163 158 DO jk = 1, jpkm1 164 zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * zrhd(:,:,jk)159 zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * rhd(:,:,jk) 165 160 END DO 166 161 IF( ln_linssh ) THEN … … 169 164 DO jj = 1,jpj 170 165 iks = mikt(ji,jj) 171 zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj)166 zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * rhd(ji,jj,iks) + riceload(ji,jj) 172 167 END DO 173 168 END DO 174 169 ELSE 175 zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1)170 zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * rhd(:,:,1) 176 171 END IF 177 172 END IF 178 173 ! 179 zarho = glob_sum( 'diaar5', area(:,:) * zbotpres(:,:) )174 zarho = glob_sum( 'diaar5', e1e2t(:,:) * zbotpres(:,:) ) 180 175 zssh_steric = - zarho / area_tot 181 176 CALL iom_put( 'sshsteric', zssh_steric ) … … 191 186 ztsn(:,:,:,:) = 0._wp ! ztsn(:,:,1,jp_tem/sal) is used here as 2D Workspace for temperature & salinity 192 187 DO_3D_11_11( 1, jpkm1 ) 193 zztmp = area(ji,jj) * e3t(ji,jj,jk,Kmm)188 zztmp = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) 194 189 ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zztmp * ts(ji,jj,jk,jp_tem,Kmm) 195 190 ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zztmp * ts(ji,jj,jk,jp_sal,Kmm) … … 237 232 z2d(:,:) = 0._wp 238 233 DO jk = 1, jpkm1 239 z2d(:,:) = z2d(:,:) + area(:,:) * e3t(:,:,jk,Kmm) * ztpot(:,:,jk)234 z2d(:,:) = z2d(:,:) + e1e2t(:,:) * e3t(:,:,jk,Kmm) * ztpot(:,:,jk) 240 235 END DO 241 236 ztemp = glob_sum( 'diaar5', z2d(:,:) ) … … 244 239 ! 245 240 IF( iom_use( 'ssttot' ) ) THEN ! Output potential temperature in case we use TEOS-10 246 zsst = glob_sum( 'diaar5', area(:,:) * ztpot(:,:,1) )241 zsst = glob_sum( 'diaar5', e1e2t(:,:) * ztpot(:,:,1) ) 247 242 CALL iom_put( 'ssttot', zsst / area_tot ) 248 243 ENDIF … … 259 254 ELSE 260 255 IF( iom_use('ssttot') ) THEN ! Output sst in case we use EOS-80 261 zsst = glob_sum( 'diaar5', area(:,:) * ts(:,:,1,jp_tem,Kmm) )256 zsst = glob_sum( 'diaar5', e1e2t(:,:) * ts(:,:,1,jp_tem,Kmm) ) 262 257 CALL iom_put('ssttot', zsst / area_tot ) 263 258 ENDIF … … 294 289 IF( l_ar5 ) THEN 295 290 DEALLOCATE( zarea_ssh , zbotpres, z2d ) 296 DEALLOCATE( zrhd , zrhop )297 291 DEALLOCATE( ztsn ) 298 292 ENDIF … … 368 362 IF( iom_use( 'voltot' ) .OR. iom_use( 'sshtot' ) .OR. iom_use( 'sshdyn' ) .OR. & 369 363 & iom_use( 'masstot' ) .OR. iom_use( 'temptot' ) .OR. iom_use( 'saltot' ) .OR. & 370 & iom_use( 'botpres' ) .OR. iom_use( 'sshthster' ) .OR. iom_use( 'sshsteric' ) ) L_ar5 = .TRUE. 364 & iom_use( 'botpres' ) .OR. iom_use( 'sshthster' ) .OR. iom_use( 'sshsteric' ) .OR. & 365 & iom_use( 'rhop' ) ) L_ar5 = .TRUE. 371 366 372 367 IF( l_ar5 ) THEN … … 375 370 IF( dia_ar5_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' ) 376 371 377 area(:,:) = e1e2t(:,:) 378 area_tot = glob_sum( 'diaar5', area(:,:) ) 372 area_tot = glob_sum( 'diaar5', e1e2t(:,:) ) 379 373 380 374 ALLOCATE( zvol0(jpi,jpj) ) … … 383 377 DO_3D_11_11( 1, jpkm1 ) 384 378 idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk) 385 zvol0 (ji,jj) = zvol0 (ji,jj) + idep * area(ji,jj)379 zvol0 (ji,jj) = zvol0 (ji,jj) + idep * e1e2t(ji,jj) 386 380 thick0(ji,jj) = thick0(ji,jj) + idep 387 381 END_3D -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DIA/diamlr.F90
r12377 r13159 84 84 INTEGER :: itide ! Number of available tidal components 85 85 REAL(wp) :: ztide_phase ! Tidal-constituent phase at adatrj=0 86 CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: ctide_selected = ' 86 CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: ctide_selected = 'n/a ' 87 87 TYPE(tide_harmonic), DIMENSION(:), POINTER :: stideconst 88 88 … … 145 145 ! Retrieve information (frequency, phase, nodal correction) about all 146 146 ! available tidal constituents for placeholder substitution below 147 ctide_selected(1:34) = (/ 'Mf', 'Mm', 'Ssa', 'Mtm', 'Msf', & 148 & 'Msqm', 'Sa', 'K1', 'O1', 'P1', & 149 & 'Q1', 'J1', 'S1', 'M2', 'S2', 'N2', & 150 & 'K2', 'nu2', 'mu2', '2N2', 'L2', & 151 & 'T2', 'eps2', 'lam2', 'R2', 'M3', & 152 & 'MKS2', 'MN4', 'MS4', 'M4', 'N4', & 153 & 'S4', 'M6', 'M8' /) 147 ! Warning: we must use the same character length in an array constructor (at least for gcc compiler) 148 ctide_selected(1:34) = (/ 'Mf ', 'Mm ', 'Ssa ', 'Mtm ', 'Msf ', & 149 & 'Msqm', 'Sa ', 'K1 ', 'O1 ', 'P1 ', & 150 & 'Q1 ', 'J1 ', 'S1 ', 'M2 ', 'S2 ', 'N2 ', & 151 & 'K2 ', 'nu2 ', 'mu2 ', '2N2 ', 'L2 ', & 152 & 'T2 ', 'eps2', 'lam2', 'R2 ', 'M3 ', & 153 & 'MKS2', 'MN4 ', 'MS4 ', 'M4 ', 'N4 ', & 154 & 'S4 ', 'M6 ', 'M8 ' /) 154 155 CALL tide_init_harmonics(ctide_selected, stideconst) 155 156 itide = size(stideconst) -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DIA/diawri.F90
r12493 r13159 171 171 CALL iom_put( "sbs", z2d ) ! bottom salinity 172 172 ENDIF 173 174 CALL iom_put( "rhop", rhop(:,:,:) ) ! 3D potential density (sigma0) 173 175 174 176 IF ( iom_use("taubot") ) THEN ! bottom stress … … 924 926 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created ' 925 927 IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc' 926 927 #if defined key_si3 928 CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 929 #else 930 CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 931 #endif 932 928 ! 929 CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 930 ! 933 931 CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) ) ! now temperature 934 932 CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) ) ! now salinity … … 943 941 CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep ) ! now k-velocity 944 942 CALL iom_rstput( 0, 0, inum, 'ht' , ht ) ! now water column height 945 943 ! 946 944 IF ( ln_isf ) THEN 947 945 IF (ln_isfcav_mlt) THEN … … 949 947 CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav ) ! now k-velocity 950 948 CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav ) ! now k-velocity 951 CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav, 8)) ! now k-velocity952 CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav, 8)) ! now k-velocity953 CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav, 8), ktype = jp_i1 )949 CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) ) ! now k-velocity 950 CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) ) ! now k-velocity 951 CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 ) 954 952 END IF 955 953 IF (ln_isfpar_mlt) THEN 956 CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par, 8)) ! now k-velocity954 CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) ) ! now k-velocity 957 955 CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par ) ! now k-velocity 958 956 CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par ) ! now k-velocity 959 957 CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par ) ! now k-velocity 960 CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par, 8)) ! now k-velocity961 CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par, 8)) ! now k-velocity962 CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par, 8), ktype = jp_i1 )958 CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) ) ! now k-velocity 959 CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) ) ! now k-velocity 960 CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 ) 963 961 END IF 964 962 END IF 965 963 ! 966 964 IF( ALLOCATED(ahtu) ) THEN 967 965 CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point … … 993 991 CALL iom_rstput ( 0, 0, inum, "qz1_abl", tq_abl(:,:,2,nt_a,2) ) ! now first level humidity 994 992 ENDIF 995 993 ! 994 CALL iom_close( inum ) 995 ! 996 996 #if defined key_si3 997 997 IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid 998 CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' ) 998 999 CALL ice_wri_state( inum ) 999 ENDIF 1000 CALL iom_close( inum ) 1001 ENDIF 1002 ! 1000 1003 #endif 1001 !1002 CALL iom_close( inum )1003 !1004 1004 END SUBROUTINE dia_wri_state 1005 1005 -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DOM/dom_oce.F90
r12489 r13159 17 17 !!---------------------------------------------------------------------- 18 18 !! Agrif_Root : dummy function used when lk_agrif=F 19 !! Agrif_Fixed : dummy function used when lk_agrif=F 19 20 !! Agrif_CFixed : dummy function used when lk_agrif=F 20 21 !! dom_oce_alloc : dynamical allocation of dom_oce arrays … … 233 234 END FUNCTION Agrif_Root 234 235 236 INTEGER FUNCTION Agrif_Fixed() 237 Agrif_Fixed = 0 238 END FUNCTION Agrif_Fixed 239 235 240 CHARACTER(len=3) FUNCTION Agrif_CFixed() 236 241 Agrif_CFixed = '0' -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DOM/dommsk.F90
r12377 r13159 259 259 ENDIF 260 260 END DO 261 #if defined key_agrif262 IF( .NOT. AGRIF_Root() ) THEN263 IF ((nbondi == 1).OR.(nbondi == 2)) fmask(nlci-1 , : ,jk) = 0.e0 ! east264 IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1 , : ,jk) = 0.e0 ! west265 IF ((nbondj == 1).OR.(nbondj == 2)) fmask(: ,nlcj-1 ,jk) = 0.e0 ! north266 IF ((nbondj == -1).OR.(nbondj == 2)) fmask(: ,1 ,jk) = 0.e0 ! south267 ENDIF268 #endif269 261 END DO 270 262 ! -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DOM/domvvl.F90
r12489 r13159 903 903 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 904 904 905 DO ji = 1, jpi 906 DO jj = 1, jpj 907 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 908 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 909 ENDIF 910 END DO 911 END DO 905 DO_2D_11_11 906 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 907 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 908 ENDIF 909 END_2D 912 910 ! 913 911 ELSE -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DOM/istate.F90
r12489 r13159 24 24 USE dom_oce ! ocean space and time domain 25 25 USE daymod ! calendar 26 USE divhor ! horizontal divergence (div_hor routine)27 26 USE dtatsd ! data temperature and salinity (dta_tsd routine) 28 27 USE dtauvd ! data: U & V current (dta_uvd routine) … … 121 120 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 122 121 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 123 hdiv(:,:,jpk) = 0._wp ! bottom divergence set one for 0 to zero at jpk level124 CALL div_hor( 0, Kbb, Kmm ) ! compute interior hdiv value125 !!gm hdiv(:,:,:) = 0._wp126 122 127 123 !!gm POTENTIAL BUG : -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DYN/divhor.F90
r12377 r13159 84 84 END_3D 85 85 ! 86 #if defined key_agrif87 IF( .NOT. Agrif_Root() ) THEN88 IF( nbondi == -1 .OR. nbondi == 2 ) hdiv( 2 , : ,:) = 0._wp ! west89 IF( nbondi == 1 .OR. nbondi == 2 ) hdiv( nlci-1, : ,:) = 0._wp ! east90 IF( nbondj == -1 .OR. nbondj == 2 ) hdiv( : , 2 ,:) = 0._wp ! south91 IF( nbondj == 1 .OR. nbondj == 2 ) hdiv( : ,nlcj-1,:) = 0._wp ! north92 ENDIF93 #endif94 !95 86 IF( ln_rnf ) CALL sbc_rnf_div( hdiv, Kmm ) !== runoffs ==! (update hdiv field) 96 87 ! -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DYN/dynldf_lap_blp.F90
r12377 r13159 74 74 DO_2D_01_01 75 75 ! ! ahm * e3 * curl (computed from 1 to jpim1/jpjm1) 76 !!gm open question here : e3f at before or now ? probably now... 77 !!gm note that ahmf has already been multiplied by fmask 78 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & 76 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & ! ahmf already * by fmask 79 77 & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) & 80 78 & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) 81 79 ! ! ahm * div (computed from 2 to jpi/jpj) 82 !!gm note that ahmt has already been multiplied by tmask 83 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & 80 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & ! ahmt already * by tmask 84 81 & * ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) & 85 82 & + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) ) … … 87 84 ! 88 85 DO_2D_00_00 89 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * ( &86 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * ( & ! * by umask is mandatory for dyn_ldf_blp use 90 87 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 91 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) )88 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) 92 89 ! 93 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * ( &90 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * ( & ! * by vmask is mandatory for dyn_ldf_blp use 94 91 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm) & 95 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) )92 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) 96 93 END_2D 97 94 ! ! =============== -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DYN/dynvor.F90
r12377 r13159 810 810 DO_3D_10_10( 1, jpk ) 811 811 IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 812 & + tmask(ji,jj ,jk) + tmask(ji+1,jj +1,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp812 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp 813 813 END_3D 814 814 ! -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DYN/sshwzv.F90
r12489 r13159 202 202 #if defined key_agrif 203 203 IF( .NOT. AGRIF_Root() ) THEN 204 IF ((nbondi == 1).OR.(nbondi == 2)) pww(nlci-1 , : ,:) = 0.e0 ! east 205 IF ((nbondi == -1).OR.(nbondi == 2)) pww(2 , : ,:) = 0.e0 ! west 206 IF ((nbondj == 1).OR.(nbondj == 2)) pww(: ,nlcj-1 ,:) = 0.e0 ! north 207 IF ((nbondj == -1).OR.(nbondj == 2)) pww(: ,2 ,:) = 0.e0 ! south 204 ! Mask vertical velocity at first/last columns/row 205 ! inside computational domain (cosmetic) 206 ! --- West --- ! 207 DO ji = mi0(2), mi1(2) 208 DO jj = 1, jpj 209 pww(ji,jj,:) = 0._wp 210 ENDDO 211 ENDDO 212 ! 213 ! --- East --- ! 214 DO ji = mi0(jpiglo-1), mi1(jpiglo-1) 215 DO jj = 1, jpj 216 pww(ji,jj,:) = 0._wp 217 ENDDO 218 ENDDO 219 ! 220 ! --- South --- ! 221 DO jj = mj0(2), mj1(2) 222 DO ji = 1, jpi 223 pww(ji,jj,:) = 0._wp 224 ENDDO 225 ENDDO 226 ! 227 ! --- North --- ! 228 DO jj = mj0(jpjglo-1), mj1(jpjglo-1) 229 DO ji = 1, jpi 230 pww(ji,jj,:) = 0._wp 231 ENDDO 232 ENDDO 208 233 ENDIF 209 234 #endif -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/FLO/floblk.F90
r12489 r13159 175 175 zgidfl(jfl) = float(iioutfl(jfl) - iiinfl(jfl)) 176 176 IF( zufl(jfl)*zuoutfl <= 0. ) THEN 177 ztxfl(jfl) = 1.E99177 ztxfl(jfl) = HUGE(1._wp) 178 178 ELSE 179 179 IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN … … 191 191 zgjdfl(jfl) = float(ijoutfl(jfl)-ijinfl(jfl)) 192 192 IF( zvfl(jfl)*zvoutfl <= 0. ) THEN 193 ztyfl(jfl) = 1.E99193 ztyfl(jfl) = HUGE(1._wp) 194 194 ELSE 195 195 IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN … … 208 208 zgkdfl(jfl) = float(ikoutfl(jfl) - ikinfl(jfl)) 209 209 IF( zwfl(jfl)*zwoutfl <= 0. ) THEN 210 ztzfl(jfl) = 1.E99210 ztzfl(jfl) = HUGE(1._wp) 211 211 ELSE 212 212 IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/ICB/icbrst.F90
r12472 r13159 188 188 ! 189 189 INTEGER :: jn ! dummy loop index 190 INTEGER :: idg ! number of digits 190 191 INTEGER :: ix_dim, iy_dim, ik_dim, in_dim 191 192 CHARACTER(len=256) :: cl_path 192 193 CHARACTER(len=256) :: cl_filename 193 CHARACTER(len=256) :: cl_kt 194 CHARACTER(len=8 ) :: cl_kt 195 CHARACTER(LEN=12 ) :: clfmt ! writing format 194 196 TYPE(iceberg), POINTER :: this 195 197 TYPE(point) , POINTER :: pt … … 211 213 ! file name 212 214 WRITE(cl_kt, '(i8.8)') kt 213 cl_filename = TRIM(cexper)//"_"// TRIM(ADJUSTL(cl_kt))//"_"//TRIM(cn_icbrst_out)215 cl_filename = TRIM(cexper)//"_"//cl_kt//"_"//TRIM(cn_icbrst_out) 214 216 IF( lk_mpp ) THEN 215 WRITE(cl_filename,'(A,"_",I4.4,".nc")') TRIM(cl_filename), narea-1 217 idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9 218 WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg ! '(a,a,ix.x,a)' 219 WRITE(cl_filename, clfmt) TRIM(cl_filename), '_', narea-1, '.nc' 216 220 ELSE 217 WRITE(cl_filename,'( A,".nc")') TRIM(cl_filename)221 WRITE(cl_filename,'(a,a)') TRIM(cl_filename), '.nc' 218 222 ENDIF 219 223 -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/ICB/icbtrj.F90
r12489 r13159 62 62 ! 63 63 INTEGER :: iret, iyear, imonth, iday 64 INTEGER :: idg ! number of digits 64 65 REAL(wp) :: zfjulday, zsec 65 66 CHARACTER(len=80) :: cl_filename 66 CHARACTER(LEN=20) :: cldate_ini, cldate_end 67 CHARACTER(LEN=12) :: clfmt ! writing format 68 CHARACTER(LEN=8 ) :: cldate_ini, cldate_end 67 69 TYPE(iceberg), POINTER :: this 68 70 TYPE(point) , POINTER :: pt … … 80 82 81 83 ! define trajectory output name 82 IF ( lk_mpp ) THEN ; WRITE(cl_filename,'("trajectory_icebergs_",A,"-",A,"_",I4.4,".nc")') & 83 & TRIM(ADJUSTL(cldate_ini)), TRIM(ADJUSTL(cldate_end)), narea-1 84 ELSE ; WRITE(cl_filename,'("trajectory_icebergs_",A,"-",A ,".nc")') & 85 & TRIM(ADJUSTL(cldate_ini)), TRIM(ADJUSTL(cldate_end)) 84 cl_filename = 'trajectory_icebergs_'//cldate_ini//'-'//cldate_end 85 IF ( lk_mpp ) THEN 86 idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9 87 WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg ! '(a,a,ix.x,a)' 88 WRITE(cl_filename, clfmt) TRIM(cl_filename), '_', narea-1, '.nc' 89 ELSE 90 WRITE(cl_filename,'(a,a)') TRIM(cl_filename), '.nc' 86 91 ENDIF 87 92 IF( lwp .AND. nn_verbose_level >= 0 ) WRITE(numout,'(2a)') 'icebergs, icb_trj_init: creating ',TRIM(cl_filename) -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/IOM/in_out_manager.F90
r12377 r13159 100 100 !!---------------------------------------------------------------------- 101 101 TYPE :: sn_ctl !: structure for control over output selection 102 LOGICAL :: l_glochk = .FALSE. !: range sanity checks are local (F) or global (T)103 ! Use global setting for debugging only;104 ! local breaches will still be reported105 ! and stop the code in most cases.106 LOGICAL :: l_allon = .FALSE. !: overall control; activate all following output options107 LOGICAL :: l_config = .FALSE. !: activate/deactivate finer control108 ! Note if l_config is True then sn_cfctl%l_allon is ignored.109 ! Otherwise setting sn_cfctl%l_allon T/F is equivalent to110 ! setting all the following logicals in this structure T/F111 ! and disabling subsetting of processors112 102 LOGICAL :: l_runstat = .FALSE. !: Produce/do not produce run.stat file (T/F) 113 103 LOGICAL :: l_trcstat = .FALSE. !: Produce/do not produce tracer.stat file (T/F) … … 169 159 INTEGER :: no_print = 0 !: optional argument of fld_fill (if present, suppress some control print) 170 160 INTEGER :: nstop = 0 !: error flag (=number of reason for a premature stop run) 161 !$AGRIF_DO_NOT_TREAT 162 INTEGER :: ngrdstop = -1 !: grid number having nstop > 1 163 !$AGRIF_END_DO_NOT_TREAT 171 164 INTEGER :: nwarn = 0 !: warning flag (=number of warning found during the run) 172 165 CHARACTER(lc) :: ctmp1, ctmp2, ctmp3 !: temporary characters 1 to 3 -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/IOM/iom.F90
r12799 r13159 111 111 CHARACTER(len=lc) :: clname 112 112 INTEGER :: irefyear, irefmonth, irefday 113 INTEGER :: ji , jkmin113 INTEGER :: ji 114 114 LOGICAL :: llrst_context ! is context related to restart 115 115 ! … … 220 220 221 221 ! Add vertical grid bounds 222 jkmin = MIN(2,jpk) ! in case jpk=1 (i.e. sas2D) 223 zt_bnds(2,: ) = gdept_1d(:) 224 zt_bnds(1,jkmin:jpk) = gdept_1d(1:jpkm1) 225 zt_bnds(1,1 ) = gdept_1d(1) - e3w_1d(1) 226 zw_bnds(1,: ) = gdepw_1d(:) 227 zw_bnds(2,1:jpkm1 ) = gdepw_1d(jkmin:jpk) 228 zw_bnds(2,jpk: ) = gdepw_1d(jpk) + e3t_1d(jpk) 222 zt_bnds(2,: ) = gdept_1d(:) 223 zt_bnds(1,2:jpk ) = gdept_1d(1:jpkm1) 224 zt_bnds(1,1 ) = gdept_1d(1) - e3w_1d(1) 225 zw_bnds(1,: ) = gdepw_1d(:) 226 zw_bnds(2,1:jpkm1) = gdepw_1d(2:jpk) 227 zw_bnds(2,jpk: ) = gdepw_1d(jpk) + e3t_1d(jpk) 229 228 CALL iom_set_axis_attr( "deptht", bounds=zw_bnds ) 230 229 CALL iom_set_axis_attr( "depthu", bounds=zw_bnds ) … … 665 664 666 665 667 SUBROUTINE iom_open( cdname, kiomid, ldwrt, kdom, ldstop, ldiof, kdlev )666 SUBROUTINE iom_open( cdname, kiomid, ldwrt, kdom, ldstop, ldiof, kdlev, cdcomp ) 668 667 !!--------------------------------------------------------------------- 669 668 !! *** SUBROUTINE iom_open *** … … 678 677 LOGICAL , INTENT(in ), OPTIONAL :: ldiof ! Interp On the Fly, needed for AGRIF (default = .FALSE.) 679 678 INTEGER , INTENT(in ), OPTIONAL :: kdlev ! number of vertical levels 679 CHARACTER(len=3), INTENT(in ), OPTIONAL :: cdcomp ! name of component calling iom_nf90_open 680 680 ! 681 681 CHARACTER(LEN=256) :: clname ! the name of the file based on cdname [[+clcpu]+clcpu] … … 823 823 ENDIF 824 824 IF( istop == nstop ) THEN ! no error within this routine 825 CALL iom_nf90_open( clname, kiomid, llwrt, llok, idompar, kdlev = kdlev )825 CALL iom_nf90_open( clname, kiomid, llwrt, llok, idompar, kdlev = kdlev, cdcomp = cdcomp ) 826 826 ENDIF 827 827 ! -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/IOM/iom_def.F90
r12377 r13159 33 33 INTEGER, PARAMETER, PUBLIC :: jpmax_vars = 1200 !: maximum number of variables in one file 34 34 INTEGER, PARAMETER, PUBLIC :: jpmax_dims = 4 !: maximum number of dimensions for one variable 35 INTEGER, PARAMETER, PUBLIC :: jpmax_digits = 5!: maximum number of digits for the cpu number in the file name35 INTEGER, PARAMETER, PUBLIC :: jpmax_digits = 9 !: maximum number of digits for the cpu number in the file name 36 36 37 37 … … 50 50 TYPE, PUBLIC :: file_descriptor 51 51 CHARACTER(LEN=240) :: name !: name of the file 52 CHARACTER(LEN=3 ) :: comp !: name of component opening the file ('OCE', 'ICE'...) 52 53 INTEGER :: nfid !: identifier of the file (0 if closed) 53 54 !: jpioipsl option has been removed) … … 64 65 REAL(kind=wp), DIMENSION(jpmax_vars) :: scf !: scale_factor of the variables 65 66 REAL(kind=wp), DIMENSION(jpmax_vars) :: ofs !: add_offset of the variables 66 INTEGER :: nlev ! number of vertical levels67 67 END TYPE file_descriptor 68 68 TYPE(file_descriptor), DIMENSION(jpmax_files), PUBLIC :: iom_file !: array containing the info for all opened files -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/IOM/iom_nf90.F90
r12377 r13159 19 19 !!---------------------------------------------------------------------- 20 20 USE dom_oce ! ocean space and time domain 21 USE sbc_oce, ONLY: jpka,ght_abl ! abl vertical level number and height21 USE sbc_oce, ONLY: ght_abl ! abl vertical level number and height 22 22 USE lbclnk ! lateal boundary condition / mpp exchanges 23 23 USE iom_def ! iom variables definitions … … 46 46 CONTAINS 47 47 48 SUBROUTINE iom_nf90_open( cdname, kiomid, ldwrt, ldok, kdompar, kdlev )48 SUBROUTINE iom_nf90_open( cdname, kiomid, ldwrt, ldok, kdompar, kdlev, cdcomp ) 49 49 !!--------------------------------------------------------------------- 50 50 !! *** SUBROUTINE iom_open *** … … 58 58 INTEGER, DIMENSION(2,5), INTENT(in ), OPTIONAL :: kdompar ! domain parameters: 59 59 INTEGER , INTENT(in ), OPTIONAL :: kdlev ! size of the ice/abl third dimension 60 CHARACTER(len=3) , INTENT(in ), OPTIONAL :: cdcomp ! name of component calling iom_nf90_open 60 61 61 62 CHARACTER(LEN=256) :: clinfo ! info character 62 63 CHARACTER(LEN=256) :: cltmp ! temporary character 64 CHARACTER(LEN=12 ) :: clfmt ! writing format 65 CHARACTER(LEN=3 ) :: clcomp ! name of component calling iom_nf90_open 66 INTEGER :: idg ! number of digits 63 67 INTEGER :: iln ! lengths of character 64 68 INTEGER :: istop ! temporary storage of nstop … … 70 74 INTEGER :: ihdf5 ! local variable for retrieval of value for NF90_HDF5 71 75 LOGICAL :: llclobber ! local definition of ln_clobber 72 INTEGER :: ilevels ! vertical levels73 76 !--------------------------------------------------------------------- 74 77 ! … … 77 80 ! 78 81 ! !number of vertical levels 79 IF( PRESENT(kdlev) ) THEN ; ilevels = kdlev ! use input value (useful for sea-ice and abl) 80 ELSE ; ilevels = jpk ! by default jpk 82 IF( PRESENT(cdcomp) ) THEN 83 IF( .NOT. PRESENT(kdlev) ) CALL ctl_stop( 'iom_nf90_open: cdcomp and kdlev must both be present' ) 84 clcomp = cdcomp ! use input value 85 ELSE 86 clcomp = 'OCE' ! by default 81 87 ENDIF 82 88 ! … … 105 111 IF( ldwrt ) THEN !* the file should be open in write mode so we create it... 106 112 IF( jpnij > 1 ) THEN 107 WRITE(cltmp,'(a,a,i4.4,a)') cdname(1:iln-1), '_', narea-1, '.nc' 113 idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9 114 WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg ! '(a,a,ix.x,a)' 115 WRITE(cltmp,clfmt) cdname(1:iln-1), '_', narea-1, '.nc' 108 116 cdname = TRIM(cltmp) 109 117 ENDIF … … 125 133 CALL iom_nf90_check(NF90_SET_FILL( if90id, NF90_NOFILL, idmy ), clinfo) 126 134 ! define dimensions 127 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'x', kdompar(1,1), idmy ), clinfo) 128 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'y', kdompar(2,1), idmy ), clinfo) 129 IF( PRESENT(kdlev) ) THEN 130 IF( kdlev == jpka ) THEN 131 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'nav_lev', kdlev, idmy ), clinfo) 132 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'time_counter', NF90_UNLIMITED, idmy ), clinfo) 133 ELSE 134 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'nav_lev', jpk, idmy ), clinfo) 135 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'time_counter', NF90_UNLIMITED, idmy ), clinfo) 136 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'numcat', kdlev, idmy ), clinfo) 137 ENDIF 138 ELSE 139 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'nav_lev', jpk, idmy ), clinfo) 140 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'time_counter', NF90_UNLIMITED, idmy ), clinfo) 141 ENDIF 135 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'x', kdompar(1,1), idmy ), clinfo) 136 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'y', kdompar(2,1), idmy ), clinfo) 137 SELECT CASE (clcomp) 138 CASE ('OCE') ; CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'nav_lev', jpk, idmy ), clinfo) 139 CASE ('ICE') ; CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'numcat', kdlev, idmy ), clinfo) 140 CASE ('ABL') ; CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'nav_lev', kdlev, idmy ), clinfo) 141 CASE ('SED') ; CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'numsed', kdlev, idmy ), clinfo) 142 CASE DEFAULT ; CALL ctl_stop( 'iom_nf90_open unknown component type' ) 143 END SELECT 144 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'time_counter', NF90_UNLIMITED, idmy ), clinfo) 142 145 ! global attributes 143 146 CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_number_total' , jpnij ), clinfo) … … 165 168 ENDDO 166 169 iom_file(kiomid)%name = TRIM(cdname) 170 iom_file(kiomid)%comp = clcomp 167 171 iom_file(kiomid)%nfid = if90id 168 172 iom_file(kiomid)%nvars = 0 169 173 iom_file(kiomid)%irec = -1 ! useless for NetCDF files, used to know if the file is in define mode 170 iom_file(kiomid)%nlev = ilevels171 174 CALL iom_nf90_check(NF90_Inquire(if90id, unlimitedDimId = iom_file(kiomid)%iduld), clinfo) 172 175 IF( iom_file(kiomid)%iduld .GE. 0 ) THEN … … 529 532 INTEGER, DIMENSION(4) :: idimid ! dimensions id 530 533 CHARACTER(LEN=256) :: clinfo ! info character 531 CHARACTER(LEN= 12), DIMENSION(5) :: cltmp ! temporary character532 534 INTEGER :: if90id ! nf90 file identifier 533 INTEGER :: idmy ! dummy variable534 535 INTEGER :: itype ! variable type 535 536 INTEGER, DIMENSION(4) :: ichunksz ! NetCDF4 chunk sizes. Will be computed using … … 540 541 ! ! when appropriate (currently chunking is applied to 4d fields only) 541 542 INTEGER :: idlv ! local variable 542 INTEGER :: idim3 ! id of the third dimension543 543 !--------------------------------------------------------------------- 544 544 ! … … 554 554 ENDIF 555 555 ! define the dimension variables if it is not already done 556 ! Warning: we must use the same character length in an array constructor (at least for gcc compiler) 557 cltmp = (/ 'nav_lon ', 'nav_lat ', 'nav_lev ', 'time_counter', 'numcat ' /) 558 CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(1)), NF90_FLOAT , (/ 1, 2 /), iom_file(kiomid)%nvid(1) ), clinfo) 559 CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(2)), NF90_FLOAT , (/ 1, 2 /), iom_file(kiomid)%nvid(2) ), clinfo) 560 CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(3)), NF90_FLOAT , (/ 3 /), iom_file(kiomid)%nvid(3) ), clinfo) 561 CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(4)), NF90_DOUBLE, (/ 4 /), iom_file(kiomid)%nvid(4) ), clinfo) 556 DO jd = 1, 2 557 CALL iom_nf90_check(NF90_INQUIRE_DIMENSION(if90id,jd,iom_file(kiomid)%cn_var(jd),iom_file(kiomid)%dimsz(jd,jd)),clinfo) 558 CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(iom_file(kiomid)%cn_var(jd)), NF90_FLOAT , (/ 1, 2 /), & 559 & iom_file(kiomid)%nvid(jd) ), clinfo) 560 END DO 561 iom_file(kiomid)%dimsz(2,1) = iom_file(kiomid)%dimsz(2,2) ! second dim of first variable 562 iom_file(kiomid)%dimsz(1,2) = iom_file(kiomid)%dimsz(1,1) ! first dim of second variable 563 DO jd = 3, 4 564 CALL iom_nf90_check(NF90_INQUIRE_DIMENSION(if90id,jd,iom_file(kiomid)%cn_var(jd),iom_file(kiomid)%dimsz(1,jd)), clinfo) 565 CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(iom_file(kiomid)%cn_var(jd)), NF90_FLOAT , (/ jd /), & 566 & iom_file(kiomid)%nvid(jd) ), clinfo) 567 END DO 562 568 ! update informations structure related the dimension variable we just added... 563 569 iom_file(kiomid)%nvars = 4 564 570 iom_file(kiomid)%luld(1:4) = (/ .FALSE., .FALSE., .FALSE., .TRUE. /) 565 iom_file(kiomid)%cn_var(1:4) = cltmp(1:4)566 571 iom_file(kiomid)%ndims(1:4) = (/ 2, 2, 1, 1 /) 567 IF( NF90_INQ_DIMID( if90id, 'numcat', idmy ) == nf90_noerr ) THEN ! add a 5th variable corresponding to the 5th dimension568 CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(5)), NF90_FLOAT , (/ 5 /), iom_file(kiomid)%nvid(5) ), clinfo)569 iom_file(kiomid)%nvars = 5570 iom_file(kiomid)%luld(5) = .FALSE.571 iom_file(kiomid)%cn_var(5) = cltmp(5)572 iom_file(kiomid)%ndims(5) = 1573 ENDIF574 ! trick: defined to 0 to say that dimension variables are defined but not yet written575 iom_file(kiomid)%dimsz(1, 1) = 0576 572 IF(lwp) WRITE(numout,*) TRIM(clinfo)//' define dimension variables done' 577 573 ENDIF … … 594 590 IF( PRESENT(pv_r0d) ) THEN ; idims = 0 595 591 ELSEIF( PRESENT(pv_r1d) ) THEN 596 IF(( SIZE(pv_r1d,1) == jpk ).OR.( SIZE(pv_r1d,1) == jpka )) THEN ; idim3 = 3 597 ELSE ; idim3 = 5 598 ENDIF 599 idims = 2 ; idimid(1:idims) = (/idim3,4/) 600 ELSEIF( PRESENT(pv_r2d) ) THEN ; idims = 3 ; idimid(1:idims) = (/1,2 ,4/) 592 idims = 2 ; idimid(1:idims) = (/3,4/) 593 ELSEIF( PRESENT(pv_r2d) ) THEN ; idims = 3 ; idimid(1:idims) = (/1,2,4/) 601 594 ELSEIF( PRESENT(pv_r3d) ) THEN 602 IF(( SIZE(pv_r3d,3) == jpk ).OR.( SIZE(pv_r3d,3) == jpka )) THEN ; idim3 = 3 603 ELSE ; idim3 = 5 604 ENDIF 605 idims = 4 ; idimid(1:idims) = (/1,2,idim3,4/) 595 idims = 4 ; idimid(1:idims) = (/1,2,3,4/) 606 596 ENDIF 607 597 IF( PRESENT(ktype) ) THEN ! variable external type … … 678 668 ! ============= 679 669 ! trick: is defined to 0 => dimension variable are defined but not yet written 680 IF( iom_file(kiomid)%dimsz(1, 1) == 0 ) THEN 681 CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lon' , idmy ) , clinfo ) 682 CALL iom_nf90_check( NF90_PUT_VAR ( if90id, idmy, glamt(ix1:ix2, iy1:iy2) ), clinfo ) 683 CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lat' , idmy ) , clinfo ) 684 CALL iom_nf90_check( NF90_PUT_VAR ( if90id, idmy, gphit(ix1:ix2, iy1:iy2) ), clinfo ) 685 CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lev' , idmy ), clinfo ) 686 IF (iom_file(kiomid)%nlev == jpka) THEN ; CALL iom_nf90_check( NF90_PUT_VAR ( if90id, idmy, ght_abl), clinfo ) 687 ELSE ; CALL iom_nf90_check( NF90_PUT_VAR ( if90id, idmy, gdept_1d), clinfo ) 688 ENDIF 689 IF( NF90_INQ_VARID( if90id, 'numcat', idmy ) == nf90_noerr ) THEN 690 CALL iom_nf90_check( NF90_PUT_VAR ( if90id, idmy, (/ (idlv, idlv = 1,iom_file(kiomid)%nlev) /)), clinfo ) 691 ENDIF 692 ! +++ WRONG VALUE: to be improved but not really useful... 693 CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'time_counter', idmy ), clinfo ) 694 CALL iom_nf90_check( NF90_PUT_VAR( if90id, idmy, kt ), clinfo ) 695 ! update the values of the variables dimensions size 696 CALL iom_nf90_check( NF90_INQUIRE_DIMENSION( if90id, 1, len = iom_file(kiomid)%dimsz(1,1) ), clinfo ) 697 CALL iom_nf90_check( NF90_INQUIRE_DIMENSION( if90id, 2, len = iom_file(kiomid)%dimsz(2,1) ), clinfo ) 698 iom_file(kiomid)%dimsz(1:2, 2) = iom_file(kiomid)%dimsz(1:2, 1) 699 CALL iom_nf90_check( NF90_INQUIRE_DIMENSION( if90id, 3, len = iom_file(kiomid)%dimsz(1,3) ), clinfo ) 700 iom_file(kiomid)%dimsz(1 , 4) = 1 ! unlimited dimension 670 IF( iom_file(kiomid)%dimsz(1, 4) == 0 ) THEN ! time_counter = 0 671 CALL iom_nf90_check( NF90_PUT_VAR( if90id, 1, glamt(ix1:ix2, iy1:iy2) ), clinfo ) 672 CALL iom_nf90_check( NF90_PUT_VAR( if90id, 2, gphit(ix1:ix2, iy1:iy2) ), clinfo ) 673 SELECT CASE (iom_file(kiomid)%comp) 674 CASE ('OCE') 675 CALL iom_nf90_check( NF90_PUT_VAR( if90id, 3, gdept_1d ), clinfo ) 676 CASE ('ABL') 677 CALL iom_nf90_check( NF90_PUT_VAR( if90id, 3, ght_abl ), clinfo ) 678 CASE DEFAULT 679 CALL iom_nf90_check( NF90_PUT_VAR( if90id, 3, (/ (idlv, idlv = 1,iom_file(kiomid)%dimsz(1,3)) /) ), clinfo ) 680 END SELECT 681 ! "wrong" value: to be improved but not really useful... 682 CALL iom_nf90_check( NF90_PUT_VAR( if90id, 4, kt ), clinfo ) 683 ! update the size of the variable corresponding to the unlimited dimension 684 iom_file(kiomid)%dimsz(1, 4) = 1 ! so we don't enter this IF case any more... 701 685 IF(lwp) WRITE(numout,*) TRIM(clinfo)//' write dimension variables done' 702 686 ENDIF -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/ISF/isfdiags.F90
r12340 r13159 88 88 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl, pfrac ! thickness of the tbl and fraction of last cell affected by the tbl 89 89 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pvar2d ! 2d var to map in 3d 90 CHARACTER(LEN= 256), INTENT(in) :: cdvar90 CHARACTER(LEN=*), INTENT(in) :: cdvar 91 91 !!--------------------------------------------------------------------- 92 92 INTEGER :: ji, jj, jk ! loop indices -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/LBC/lib_mpp.F90
r12512 r13159 1112 1112 CHARACTER(len=*), INTENT(in ), OPTIONAL :: cd2, cd3, cd4, cd5 1113 1113 CHARACTER(len=*), INTENT(in ), OPTIONAL :: cd6, cd7, cd8, cd9, cd10 1114 ! 1115 CHARACTER(LEN=8) :: clfmt ! writing format 1116 INTEGER :: inum 1114 1117 !!---------------------------------------------------------------------- 1115 1118 ! 1116 1119 nstop = nstop + 1 1117 1120 ! 1118 ! force to open ocean.output file if not already opened 1119 IF( numout == 6 ) CALL ctl_opn( numout, 'ocean.output', 'APPEND', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. ) 1121 IF( cd1 == 'STOP' .AND. narea /= 1 ) THEN ! Immediate stop: add an arror message in 'ocean.output' file 1122 CALL ctl_opn( inum, 'ocean.output', 'APPEND', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. ) 1123 WRITE(inum,*) 1124 WRITE(inum,*) ' ==>>> Look for "E R R O R" messages in all existing *ocean.output* files' 1125 CLOSE(inum) 1126 ENDIF 1127 IF( numout == 6 ) THEN ! force to open ocean.output file if not already opened 1128 CALL ctl_opn( numout, 'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE., narea ) 1129 ENDIF 1120 1130 ! 1121 1131 WRITE(numout,*) … … 1145 1155 WRITE(numout,*) 'huge E-R-R-O-R : immediate stop' 1146 1156 WRITE(numout,*) 1157 CALL FLUSH(numout) 1158 CALL SLEEP(60) ! make sure that all output and abort files are written by all cores. 60s should be enough... 1147 1159 CALL mppstop( ld_abort = .true. ) 1148 1160 ENDIF … … 1207 1219 ! 1208 1220 CHARACTER(len=80) :: clfile 1221 CHARACTER(LEN=10) :: clfmt ! writing format 1209 1222 INTEGER :: iost 1223 INTEGER :: idg ! number of digits 1210 1224 !!---------------------------------------------------------------------- 1211 1225 ! … … 1214 1228 clfile = TRIM(cdfile) 1215 1229 IF( PRESENT( karea ) ) THEN 1216 IF( karea > 1 ) WRITE(clfile, "(a,'_',i4.4)") TRIM(clfile), karea-1 1230 IF( karea > 1 ) THEN 1231 ! Warning: jpnij is maybe not already defined when calling ctl_opn -> use mppsize instead of jpnij 1232 idg = MAX( INT(LOG10(REAL(MAX(1,mppsize-1),wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9 1233 WRITE(clfmt, "('(a,a,i', i1, '.', i1, ')')") idg, idg ! '(a,a,ix.x)' 1234 WRITE(clfile, clfmt) TRIM(clfile), '_', karea-1 1235 ENDIF 1217 1236 ENDIF 1218 1237 #if defined key_agrif -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/LBC/mpp_loc_generic.h90
r10716 r13159 32 32 REAL(wp) , INTENT( out) :: pmin ! Global minimum of ptab 33 33 INDEX_TYPE(:) ! index of minimum in global frame 34 # if defined key_mpp_mpi35 34 ! 36 35 INTEGER :: ierror, ii, idim … … 56 55 ! 57 56 kindex(1) = mig( ilocs(1) ) 58 # 57 #if defined DIM_2d || defined DIM_3d /* avoid warning when kindex has 1 element */ 59 58 kindex(2) = mjg( ilocs(2) ) 60 # 61 # 59 #endif 60 #if defined DIM_3d /* avoid warning when kindex has 2 elements */ 62 61 kindex(3) = ilocs(3) 63 # 62 #endif 64 63 ! 65 64 DEALLOCATE (ilocs) 66 65 ! 67 66 index0 = kindex(1)-1 ! 1d index starting at 0 68 # 67 #if defined DIM_2d || defined DIM_3d /* avoid warning when kindex has 1 element */ 69 68 index0 = index0 + jpiglo * (kindex(2)-1) 70 # 71 # 69 #endif 70 #if defined DIM_3d /* avoid warning when kindex has 2 elements */ 72 71 index0 = index0 + jpiglo * jpjglo * (kindex(3)-1) 73 # 72 #endif 74 73 END IF 75 74 zain(1,:) = zmin … … 77 76 ! 78 77 IF( ln_timing ) CALL tic_tac(.TRUE., ld_global = .TRUE.) 78 #if defined key_mpp_mpi 79 79 CALL MPI_ALLREDUCE( zain, zaout, 1, MPI_2DOUBLE_PRECISION, MPI_OPERATION ,MPI_COMM_OCE, ierror) 80 #else 81 zaout(:,:) = zain(:,:) 82 #endif 80 83 IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 81 84 ! 82 85 pmin = zaout(1,1) 83 86 index0 = NINT( zaout(2,1) ) 84 # 87 #if defined DIM_3d /* avoid warning when kindex has 2 elements */ 85 88 kindex(3) = index0 / (jpiglo*jpjglo) 86 89 index0 = index0 - kindex(3) * (jpiglo*jpjglo) 87 # 88 # 90 #endif 91 #if defined DIM_2d || defined DIM_3d /* avoid warning when kindex has 1 element */ 89 92 kindex(2) = index0 / jpiglo 90 93 index0 = index0 - kindex(2) * jpiglo 91 # 94 #endif 92 95 kindex(1) = index0 93 96 kindex(:) = kindex(:) + 1 ! start indices at 1 94 #else95 kindex = 0 ; pmin = 0.96 WRITE(*,*) 'ROUTINE_LOC: You should not have seen this print! error?'97 #endif98 97 99 98 END SUBROUTINE ROUTINE_LOC -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/OBS/obs_grid.F90
r10068 r13159 684 684 & fhistx1, fhistx2, fhisty1, fhisty2 685 685 REAL(wp) :: histtol 686 686 CHARACTER(LEN=26) :: clfmt ! writing format 687 INTEGER :: idg ! number of digits 688 687 689 IF (ln_grid_search_lookup) THEN 688 690 … … 709 711 710 712 IF ( ln_grid_global ) THEN 711 WRITE(cfname, FMT="(A,'_',A)") & 712 & TRIM(cn_gridsearchfile), 'global.nc' 713 WRITE(cfname, FMT="(A,'_',A)") TRIM(cn_gridsearchfile), 'global.nc' 713 714 ELSE 714 WRITE(cfname, FMT="(A,'_',I4.4,'of',I4.4,'by',I4.4,'.nc')") & 715 & TRIM(cn_gridsearchfile), nproc, jpni, jpnj 715 idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9 716 ! define the following format: "(a,a,ix.x,a,ix.x,a,ix.x,a)" 717 WRITE(clfmt, "('(a,a,i', i1, '.', i1',a,i', i1, '.', i1',a,i', i1, '.', i1',a)')") idg, idg, idg, idg, idg, idg 718 WRITE(cfname, clfmt ) TRIM(cn_gridsearchfile),'_', nproc,'of', jpni,'by', jpnj,'.nc' 716 719 ENDIF 717 720 -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/OBS/obs_write.F90
r12377 r13159 86 86 CHARACTER(LEN=40) :: clfname 87 87 CHARACTER(LEN=10) :: clfiletype 88 CHARACTER(LEN=12) :: clfmt ! writing format 89 INTEGER :: idg ! number of digits 88 90 INTEGER :: ilevel 89 91 INTEGER :: jvar … … 181 183 fbdata%caddname(1) = 'Hx' 182 184 183 WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc 185 idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9 186 WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg ! '(a,a,ix.x,a)' 187 WRITE(clfname,clfmt) TRIM(clfiletype), '_fdbk_', nproc, '.nc' 184 188 185 189 IF(lwp) THEN … … 326 330 CHARACTER(LEN=10) :: clfiletype 327 331 CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 332 CHARACTER(LEN=12) :: clfmt ! writing format 333 INTEGER :: idg ! number of digits 328 334 INTEGER :: jo 329 335 INTEGER :: ja … … 453 459 fbdata%caddname(1) = 'Hx' 454 460 455 WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc 461 idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9 462 WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg ! '(a,a,ix.x,a)' 463 WRITE(clfname,clfmt) TRIM(clfiletype), '_fdbk_', nproc, '.nc' 456 464 457 465 IF(lwp) THEN -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90
r12799 r13159 632 632 633 633 END SELECT 634 634 635 IF( iom_use('Cd_oce') ) CALL iom_put("Cd_oce", zcd_oce * tmask(:,:,1)) 636 IF( iom_use('Ce_oce') ) CALL iom_put("Ce_oce", zce_oce * tmask(:,:,1)) 637 IF( iom_use('Ch_oce') ) CALL iom_put("Ch_oce", zch_oce * tmask(:,:,1)) 638 !! LB: mainly here for debugging purpose: 639 IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ztpot-rt0) * tmask(:,:,1)) ! potential temperature at z=zt 640 IF( iom_use('q_zt') ) CALL iom_put("q_zt", zqair * tmask(:,:,1)) ! specific humidity " 641 IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (t_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu 642 IF( iom_use('q_zu') ) CALL iom_put("q_zu", q_zu * tmask(:,:,1)) ! specific humidity " 643 IF( iom_use('ssq') ) CALL iom_put("ssq", pssq * tmask(:,:,1)) ! saturation specific humidity at z=0 644 IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu 645 635 646 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 636 647 !! ptsk and pssq have been updated!!! … … 643 654 END WHERE 644 655 END IF 645 646 !! CALL iom_put( "Cd_oce", zcd_oce) ! output value of pure ocean-atm. transfer coef.647 !! CALL iom_put( "Ch_oce", zch_oce) ! output value of pure ocean-atm. transfer coef.648 649 IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN650 !! If zu == zt, then ensuring once for all that:651 t_zu(:,:) = ztpot(:,:)652 q_zu(:,:) = zqair(:,:)653 ENDIF654 655 656 656 657 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbcblk_phy.F90 … … 671 672 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 672 673 & taum(:,:), psen(:,:), zqla(:,:), & 673 & pEvap=pevp(:,:), prhoa=rhoa(:,:) 674 & pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac ) 674 675 675 676 zqla(:,:) = zqla(:,:) * tmask(:,:,1) … … 688 689 ! ... utau, vtau at U- and V_points, resp. 689 690 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 690 ! Note th e use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves691 DO_2D_ 10_10691 ! Note that coastal wind stress is not used in the code... so this extra care has no effect 692 DO_2D_00_00 692 693 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) & 693 694 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) … … 893 894 Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical 894 895 ENDIF 895 896 !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef. 897 !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef. 898 896 897 IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice) 898 IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice) 899 IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice) 900 899 901 ! local scalars ( place there for vector optimisation purposes) 900 902 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 901 903 902 904 IF( ln_blk ) THEN 903 ! ------------------------------------------------------------ ! 904 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 905 ! ------------------------------------------------------------ ! 906 ! C-grid ice dynamics : U & V-points (same as ocean) 907 DO_2D_00_00 908 putaui(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * zcd_dui(ji+1,jj) & 909 & + rhoa(ji ,jj) * zcd_dui(ji ,jj) ) & 910 & * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 911 pvtaui(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * zcd_dui(ji,jj+1) & 912 & + rhoa(ji,jj ) * zcd_dui(ji,jj ) ) & 913 & * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 905 ! ------------------------------------------------------------- ! 906 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 907 ! ------------------------------------------------------------- ! 908 zztmp1 = rn_vfac * 0.5_wp 909 DO_2D_01_01 ! at T point 910 putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * ( pwndi(ji,jj) - zztmp1 * ( puice(ji-1,jj ) + puice(ji,jj) ) ) 911 pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * ( pwndj(ji,jj) - zztmp1 * ( pvice(ji ,jj-1) + pvice(ji,jj) ) ) 912 END_2D 913 ! 914 DO_2D_00_00 ! U & V-points (same as ocean). 915 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 916 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) 917 zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) ) 918 putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj ) ) 919 pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) ) 914 920 END_2D 915 921 CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) … … 1050 1056 evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation 1051 1057 devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT 1052 zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean1058 zevap (:,:) = emp(:,:) + tprecip(:,:) ! evaporation over ocean !LB: removed rn_efac here, correct??? 1053 1059 1054 1060 ! --- evaporation minus precipitation --- ! -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk_algo_coare3p0.F90
r12377 r13159 194 194 IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P0_INIT(l_use_cs, l_use_wl) 195 195 196 l_zt_equal_zu = .FALSE. 197 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 196 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 198 197 IF( .NOT. l_zt_equal_zu ) ALLOCATE( zeta_t(jpi,jpj) ) 199 198 -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk_algo_coare3p6.F90
r12377 r13159 194 194 IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl) 195 195 196 l_zt_equal_zu = .FALSE. 197 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 196 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 198 197 IF( .NOT. l_zt_equal_zu ) ALLOCATE( zeta_t(jpi,jpj) ) 199 198 -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r12377 r13159 98 98 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 99 99 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 100 !!---------------------------------------------------------------------- 100 !!---------------------------------------------------------------------------------- 101 101 !! *** ROUTINE turb_ecmwf *** 102 102 !! … … 184 184 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 185 185 ! 186 REAL(wp), DIMENSION(jpi,jpj) :: 187 REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu 188 REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air186 REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star 187 REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu 188 REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air 189 189 REAL(wp), DIMENSION(jpi,jpj) :: Linv !: 1/L (inverse of Monin Obukhov length... 190 190 REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q … … 196 196 CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90' 197 197 !!---------------------------------------------------------------------------------- 198 199 198 IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 200 199 201 l_zt_equal_zu = .FALSE. 202 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 200 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 203 201 204 202 !! Initializations for cool skin and warm layer: -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk_algo_ncar.F90
r12377 r13159 112 112 REAL(wp), DIMENSION(jpi,jpj) :: stab ! stability test integer 113 113 !!---------------------------------------------------------------------------------- 114 ! 115 l_zt_equal_zu = .FALSE. 116 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 114 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 117 115 118 116 U_blk = MAX( 0.5_wp , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s … … 143 141 ENDIF 144 142 145 !! Initializing values at z_u with z_t values: 146 t_zu = t_zt ; q_zu = q_zt 143 !! First guess of temperature and humidity at height zu: 144 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... 145 q_zu = MAX( q_zt , 1.e-6_wp ) ! " 147 146 148 147 !! ITERATION BLOCK -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk_phy.F90
r12377 r13159 31 31 REAL(wp), PARAMETER, PUBLIC :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg] 32 32 REAL(wp), PARAMETER, PUBLIC :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622 33 REAL(wp), PARAMETER, PUBLIC :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.60833 REAL(wp), PARAMETER, PUBLIC :: rctv0 = R_vap/R_dry - 1._wp !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 34 34 REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp !: specific heat of air (only used for ice fluxes now...) 35 35 REAL(wp), PARAMETER, PUBLIC :: rCd_ice = 1.4e-3_wp !: transfer coefficient over ice … … 520 520 zCe = zz0*pqst(ji,jj)/zdq 521 521 522 CALL BULK_FORMULA ( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &523 & zCd, zCh, zCe,&524 & pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),&525 & pTau(ji,jj), zQsen, zQlat )522 CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), & 523 & zCd, zCh, zCe, & 524 & pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj), & 525 & pTau(ji,jj), zQsen, zQlat ) 526 526 527 527 zTs2 = pTs(ji,jj)*pTs(ji,jj) … … 535 535 536 536 537 SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, & 538 & pCd, pCh, pCe, & 539 & pwnd, pUb, pslp, & 540 & pTau, pQsen, pQlat, pEvap, prhoa ) 537 SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, & 538 & pCd, pCh, pCe, & 539 & pwnd, pUb, pslp, & 540 & pTau, pQsen, pQlat, & 541 & pEvap, prhoa, pfact_evap ) 542 !!---------------------------------------------------------------------------------- 543 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m) 544 REAL(wp), INTENT(in) :: pTs ! water temperature at the air-sea interface [K] 545 REAL(wp), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg] 546 REAL(wp), INTENT(in) :: pTa ! potential air temperature at z=pzu [K] 547 REAL(wp), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg] 548 REAL(wp), INTENT(in) :: pCd 549 REAL(wp), INTENT(in) :: pCh 550 REAL(wp), INTENT(in) :: pCe 551 REAL(wp), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s] 552 REAL(wp), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 553 REAL(wp), INTENT(in) :: pslp ! sea-level atmospheric pressure [Pa] 554 !! 555 REAL(wp), INTENT(out) :: pTau ! module of the wind stress [N/m^2] 556 REAL(wp), INTENT(out) :: pQsen ! [W/m^2] 557 REAL(wp), INTENT(out) :: pQlat ! [W/m^2] 558 !! 559 REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 560 REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3] 561 REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent) 562 !! 563 REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap 564 INTEGER :: jq 565 !!---------------------------------------------------------------------------------- 566 zfact_evap = 1._wp 567 IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap 568 569 !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 570 ztaa = pTa ! first guess... 571 DO jq = 1, 4 572 zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa ) !LOLO: why not "0.5*(pqs+pqa)" rather then "pqa" ??? 573 ztaa = pTa - zgamma*pzu ! Absolute temp. is slightly colder... 574 END DO 575 zrho = rho_air(ztaa, pqa, pslp) 576 zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 577 578 zUrho = pUb*MAX(zrho, 1._wp) ! rho*U10 579 580 pTau = zUrho * pCd * pwnd ! Wind stress module 581 582 zevap = zUrho * pCe * (pqa - pqs) 583 pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa) 584 pQlat = L_vap(pTs) * zevap 585 586 IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap 587 IF( PRESENT(prhoa) ) prhoa = zrho 588 589 END SUBROUTINE BULK_FORMULA_SCLR 590 591 SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, & 592 & pCd, pCh, pCe, & 593 & pwnd, pUb, pslp, & 594 & pTau, pQsen, pQlat, & 595 & pEvap, prhoa, pfact_evap ) 541 596 !!---------------------------------------------------------------------------------- 542 597 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m) … … 558 613 REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 559 614 REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3] 560 !! 561 REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap 562 INTEGER :: ji, jj, jq ! dummy loop indices 563 !!---------------------------------------------------------------------------------- 564 DO_2D_11_11 565 566 !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 567 ztaa = pTa(ji,jj) ! first guess... 568 DO jq = 1, 4 569 zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) ) 570 ztaa = pTa(ji,jj) - zgamma*pzu ! Absolute temp. is slightly colder... 571 END DO 572 zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)) 573 zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 574 575 zUrho = pUb(ji,jj)*MAX(zrho, 1._wp) ! rho*U10 576 577 pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module 578 579 zevap = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj)) 580 pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj)) 581 pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap 582 583 IF( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 615 REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent) 616 !! 617 REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap 618 INTEGER :: ji, jj 619 !!---------------------------------------------------------------------------------- 620 zfact_evap = 1._wp 621 IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap 622 623 DO_2D_11_11 624 625 CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), & 626 & pCd(ji,jj), pCh(ji,jj), pCe(ji,jj), & 627 & pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj), & 628 & pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj), & 629 & pEvap=zevap, prhoa=zrho, pfact_evap=zfact_evap ) 630 631 IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap 584 632 IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho 585 633 586 634 END_2D 587 635 END SUBROUTINE BULK_FORMULA_VCTR 588 589 590 SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &591 & pCd, pCh, pCe, &592 & pwnd, pUb, pslp, &593 & pTau, pQsen, pQlat, pEvap, prhoa )594 !!----------------------------------------------------------------------------------595 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)596 REAL(wp), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]597 REAL(wp), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]598 REAL(wp), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]599 REAL(wp), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]600 REAL(wp), INTENT(in) :: pCd601 REAL(wp), INTENT(in) :: pCh602 REAL(wp), INTENT(in) :: pCe603 REAL(wp), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]604 REAL(wp), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]605 REAL(wp), INTENT(in) :: pslp ! sea-level atmospheric pressure [Pa]606 !!607 REAL(wp), INTENT(out) :: pTau ! module of the wind stress [N/m^2]608 REAL(wp), INTENT(out) :: pQsen ! [W/m^2]609 REAL(wp), INTENT(out) :: pQlat ! [W/m^2]610 !!611 REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]612 REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]613 !!614 REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap615 INTEGER :: jq616 !!----------------------------------------------------------------------------------617 618 !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")619 ztaa = pTa ! first guess...620 DO jq = 1, 4621 zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )622 ztaa = pTa - zgamma*pzu ! Absolute temp. is slightly colder...623 END DO624 zrho = rho_air(ztaa, pqa, pslp)625 zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!626 627 zUrho = pUb*MAX(zrho, 1._wp) ! rho*U10628 629 pTau = zUrho * pCd * pwnd ! Wind stress module630 631 zevap = zUrho * pCe * (pqa - pqs)632 pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)633 pQlat = L_vap(pTs) * zevap634 635 IF( PRESENT(pEvap) ) pEvap = - zevap636 IF( PRESENT(prhoa) ) prhoa = zrho637 638 END SUBROUTINE BULK_FORMULA_SCLR639 640 641 636 642 637 -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbccpl.F90
r12489 r13159 364 364 ! 365 365 ! Vectors: change of sign at north fold ONLY if on the local grid 366 IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM(sn_rcv_tau%cldes ) == 'oce and ice') THEN ! avoid working with the atmospheric fields if they are not coupled 366 IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice' & 367 .OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-ice' ) THEN ! avoid working with the atmospheric fields if they are not coupled 368 367 369 IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' ) srcv(jpr_otx1:jpr_itz2)%nsgn = -1. 368 370 … … 1115 1117 IF( ln_dm2dc .AND. ncpl_qsr_freq /= 86400 ) & 1116 1118 & CALL ctl_stop( 'sbc_cpl_rcv: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' ) 1117 ncpl_qsr_freq = 86400 / ncpl_qsr_freq ! used by top 1119 1120 IF( ncpl_qsr_freq /= 0) ncpl_qsr_freq = 86400 / ncpl_qsr_freq ! used by top 1121 1118 1122 ENDIF 1119 1123 ! … … 1479 1483 INTEGER :: ji, jj ! dummy loop indices 1480 1484 INTEGER :: itx ! index of taux over ice 1485 REAL(wp) :: zztmp1, zztmp2 1481 1486 REAL(wp), DIMENSION(jpi,jpj) :: ztx, zty 1482 1487 !!---------------------------------------------------------------------- … … 1542 1547 p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1) ! (U,V) ==> (U,V) 1543 1548 p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1) 1544 CASE( 'F' )1545 DO_2D_00_001546 p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji ,jj-1,1) )1547 p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj ,1) )1548 END_2D1549 1549 CASE( 'T' ) 1550 1550 DO_2D_00_00 1551 p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj ,1) + frcv(jpr_itx1)%z3(ji,jj,1) ) 1552 p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) ) 1551 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 1552 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) 1553 zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) ) 1554 p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj ,1) + frcv(jpr_itx1)%z3(ji,jj,1) ) 1555 p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) ) 1553 1556 END_2D 1554 CASE( 'I' ) 1555 DO_2D_00_00 1556 p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj ,1) ) 1557 p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji ,jj+1,1) ) 1558 END_2D 1557 CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U', -1., p_tauj, 'V', -1. ) 1559 1558 END SELECT 1560 IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN1561 CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U', -1., p_tauj, 'V', -1. )1562 ENDIF1563 1559 1564 1560 ENDIF … … 1789 1785 ENDDO 1790 1786 ELSE 1791 qns_tot(:,:) =qns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)1787 zqns_tot(:,:) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1) 1792 1788 DO jl = 1, jpl 1793 zqns_tot(:,: ) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)1794 1789 zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) 1795 1790 END DO … … 1932 1927 END DO 1933 1928 ELSE 1934 qsr_tot(:,: ) =qsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)1929 zqsr_tot(:,:) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1) 1935 1930 DO jl = 1, jpl 1936 zqsr_tot(:,: ) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)1937 1931 zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1) 1938 1932 END DO -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcmod.F90
r12489 r13159 120 120 ncom_fsbc = nn_fsbc ! make nn_fsbc available for lib_mpp 121 121 #endif 122 ! !* overwrite namelist parameter using CPP key information123 #if defined key_agrif124 IF( Agrif_Root() ) THEN ! AGRIF zoom (cf r1242: possibility to run without ice in fine grid)125 IF( lk_si3 ) nn_ice = 2126 IF( lk_cice ) nn_ice = 3127 ENDIF128 !!GS: TBD129 !#else130 ! IF( lk_si3 ) nn_ice = 2131 ! IF( lk_cice ) nn_ice = 3132 #endif133 122 ! 134 123 IF(lwp) THEN !* Control print -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcwave.F90
r12377 r13159 210 210 END_3D 211 211 ! 212 #if defined key_agrif213 IF( .NOT. Agrif_Root() ) THEN214 IF( nbondi == -1 .OR. nbondi == 2 ) ze3divh( 2:nbghostcells+1,: ,:) = 0._wp ! west215 IF( nbondi == 1 .OR. nbondi == 2 ) ze3divh( nlci-nbghostcells:nlci-1,:,:) = 0._wp ! east216 IF( nbondj == -1 .OR. nbondj == 2 ) ze3divh( :,2:nbghostcells+1 ,:) = 0._wp ! south217 IF( nbondj == 1 .OR. nbondj == 2 ) ze3divh( :,nlcj-nbghostcells:nlcj-1,:) = 0._wp ! north218 ENDIF219 #endif220 !221 212 CALL lbc_lnk( 'sbcwave', ze3divh, 'T', 1. ) 222 213 ! -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/STO/stopar.F90
r12377 r13159 684 684 !! ** Purpose : read stochastic parameters from restart file 685 685 !!---------------------------------------------------------------------- 686 INTEGER :: jsto, jseed 686 INTEGER :: jsto, jseed 687 INTEGER :: idg ! number of digits 687 688 INTEGER(KIND=8) :: ziseed(4) ! RNG seeds in integer type 688 689 REAL(KIND=8) :: zrseed(4) ! RNG seeds in real type (with same bits to save in restart) 689 690 CHARACTER(LEN=9) :: clsto2d='sto2d_000' ! stochastic parameter variable name 690 691 CHARACTER(LEN=9) :: clsto3d='sto3d_000' ! stochastic parameter variable name 691 CHARACTER(LEN=10) :: clseed='seed0_0000' ! seed variable name 692 CHARACTER(LEN=15) :: clseed='seed0_0000' ! seed variable name 693 CHARACTER(LEN=6) :: clfmt ! writing format 692 694 !!---------------------------------------------------------------------- 693 695 … … 717 719 IF (ln_rstseed) THEN 718 720 ! Get saved state of the random number generator 721 idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9 722 WRITE(clfmt, "('(i', i1, '.', i1, ')')") idg, idg ! "(ix.x)" 719 723 DO jseed = 1 , 4 720 WRITE(clseed(5:5) ,'(i1.1)') jseed721 WRITE(clseed(7: 10),'(i4.4)') narea722 CALL iom_get( numstor, clseed , zrseed(jseed) )724 WRITE(clseed(5:5) ,'(i1.1)') jseed 725 WRITE(clseed(7:7+idg-1), clfmt ) narea 726 CALL iom_get( numstor, clseed(1:7+idg-1) , zrseed(jseed) ) 723 727 END DO 724 728 ziseed = TRANSFER( zrseed , ziseed) … … 742 746 INTEGER, INTENT(in) :: kt ! ocean time-step 743 747 !! 744 INTEGER :: jsto, jseed 748 INTEGER :: jsto, jseed 749 INTEGER :: idg ! number of digits 745 750 INTEGER(KIND=8) :: ziseed(4) ! RNG seeds in integer type 746 751 REAL(KIND=8) :: zrseed(4) ! RNG seeds in real type (with same bits to save in restart) … … 749 754 CHARACTER(LEN=9) :: clsto2d='sto2d_000' ! stochastic parameter variable name 750 755 CHARACTER(LEN=9) :: clsto3d='sto3d_000' ! stochastic parameter variable name 751 CHARACTER(LEN=10) :: clseed='seed0_0000' ! seed variable name 756 CHARACTER(LEN=15) :: clseed='seed0_0000' ! seed variable name 757 CHARACTER(LEN=6) :: clfmt ! writing format 752 758 !!---------------------------------------------------------------------- 753 759 … … 771 777 CALL kiss_state( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) ) 772 778 zrseed = TRANSFER( ziseed , zrseed) 779 idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9 780 WRITE(clfmt, "('(i', i1, '.', i1, ')')") idg, idg ! "(ix.x)" 773 781 DO jseed = 1 , 4 774 WRITE(clseed(5:5) ,'(i1.1)') jseed775 WRITE(clseed(7: 10),'(i4.4)') narea776 CALL iom_rstput( kt, nitrst, numstow, clseed 782 WRITE(clseed(5:5) ,'(i1.1)') jseed 783 WRITE(clseed(7:7+idg-1), clfmt ) narea 784 CALL iom_rstput( kt, nitrst, numstow, clseed(1:7+idg-1), zrseed(jseed) ) 777 785 END DO 778 786 ! 2D stochastic parameters -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/STO/storng.F90
r12377 r13159 50 50 51 51 ! Parameters to generate real random variates 52 REAL(KIND=wp), PARAMETER :: huge64=9223372036854775808.0 ! +153 52 REAL(KIND=wp), PARAMETER :: zero=0.0, half=0.5, one=1.0, two=2.0 54 53 … … 275 274 REAL(KIND=wp) :: uran 276 275 277 uran = half * ( one + REAL(kiss(),wp) / huge64)276 uran = half * ( one + REAL(kiss(),wp) / HUGE(1._wp) ) 278 277 279 278 END SUBROUTINE kiss_uniform … … 298 297 rsq = two 299 298 DO WHILE ( (rsq.GE.one).OR. (rsq.EQ.zero) ) 300 u1 = REAL(kiss(),wp) / huge64301 u2 = REAL(kiss(),wp) / huge64299 u1 = REAL(kiss(),wp) / HUGE(1._wp) 300 u2 = REAL(kiss(),wp) / HUGE(1._wp) 302 301 rsq = u1*u1 + u2*u2 303 302 ENDDO -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/TRD/trdtra.F90
r12489 r13159 82 82 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: ptra ! now tracer variable 83 83 ! 84 INTEGER :: jk ! loop indices 84 INTEGER :: jk ! loop indices 85 INTEGER :: i01 ! 0 or 1 85 86 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztrds ! 3D workspace 86 87 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwt, zws, ztrdt ! 3D workspace … … 90 91 IF( trd_tra_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_tra : unable to allocate arrays' ) 91 92 ENDIF 92 93 ! 94 i01 = COUNT( (/ PRESENT(pu) .OR. ( ktrd /= jptra_xad .AND. ktrd /= jptra_yad .AND. ktrd /= jptra_zad ) /) ) 95 ! 93 96 IF( ctype == 'TRA' .AND. ktra == jp_tem ) THEN !== Temperature trend ==! 94 97 ! 95 SELECT CASE( ktrd )98 SELECT CASE( ktrd*i01 ) 96 99 ! ! advection: transform the advective flux into a trend 97 100 CASE( jptra_xad ) ; CALL trd_tra_adv( ptrd, pu, ptra, 'X', trdtx, Kmm ) … … 112 115 IF( ctype == 'TRA' .AND. ktra == jp_sal ) THEN !== Salinity trends ==! 113 116 ! 114 SELECT CASE( ktrd )117 SELECT CASE( ktrd*i01 ) 115 118 ! ! advection: transform the advective flux into a trend 116 119 ! ! and send T & S trends to trd_tra_mng … … 163 166 IF( ctype == 'TRC' ) THEN !== passive tracer trend ==! 164 167 ! 165 SELECT CASE( ktrd )168 SELECT CASE( ktrd*i01 ) 166 169 ! ! advection: transform the advective flux into a masked trend 167 170 CASE( jptra_xad ) ; CALL trd_tra_adv( ptrd , pu , ptra, 'X', ztrds, Kmm ) -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/USR/usrdef_zgr.F90
r12377 r13159 202 202 CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! set surrounding land to zero (here jperio=0 ==>> closed) 203 203 ! 204 k_bot(:,:) = NINT( z2d(:,:) ) 204 k_bot(:,:) = NINT( z2d(:,:) ) ! =jpkm1 over the ocean point, =0 elsewhere 205 205 ! 206 206 k_top(:,:) = MIN( 1 , k_bot(:,:) ) ! = 1 over the ocean point, =0 elsewhere -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/ZDF/zdftke.F90
r12489 r13159 45 45 USE zdfdrg ! vertical physics: top/bottom drag coef. 46 46 USE zdfmxl ! vertical physics: mixed layer 47 #if defined key_si3 48 USE ice, ONLY: hm_i, h_i 49 #endif 50 #if defined key_cice 51 USE sbc_ice, ONLY: h_i 52 #endif 47 53 ! 48 54 USE in_out_manager ! I/O manager … … 64 70 INTEGER :: nn_mxl ! type of mixing length (=0/1/2/3) 65 71 REAL(wp) :: rn_mxl0 ! surface min value of mixing length (kappa*z_o=0.4*0.1 m) [m] 72 INTEGER :: nn_mxlice ! type of scaling under sea-ice 73 REAL(wp) :: rn_mxlice ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1) 66 74 INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1) 67 75 REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) … … 214 222 ! ! Surface/top/bottom boundary condition on tke 215 223 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 216 224 ! 217 225 DO_2D_00_00 218 226 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 219 227 END_2D 220 IF ( ln_isfcav ) THEN221 DO_2D_00_00222 en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1)223 END_2D224 ENDIF225 228 ! 226 229 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 249 252 zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2 & 250 253 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2 ) 251 en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) ! masked at ocean surface 254 ! (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) = 1 where ice shelves are present 255 en(ji,jj,mikt(ji,jj)) = en(ji,jj,1) * tmask(ji,jj,1) & 256 & + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) 252 257 END_2D 253 258 ENDIF … … 425 430 REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars 426 431 REAL(wp) :: zdku, zdkv, zsqen ! - - 427 REAL(wp) :: zemxl, zemlm, zemlp ! - -432 REAL(wp) :: zemxl, zemlm, zemlp, zmaxice ! - - 428 433 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxlm, zmxld ! 3D workspace 429 434 !!-------------------------------------------------------------------- … … 439 444 zmxld(:,:,:) = rmxl_min 440 445 ! 441 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 446 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 447 ! 442 448 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 449 #if ! defined key_si3 && ! defined key_cice 443 450 DO_2D_00_00 444 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1))451 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 445 452 END_2D 446 ELSE 453 #else 454 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 455 ! 456 CASE( 0 ) ! No scaling under sea-ice 457 DO_2D_00_00 458 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 459 END_2D 460 ! 461 CASE( 1 ) ! scaling with constant sea-ice thickness 462 DO_2D_00_00 463 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 464 END_2D 465 ! 466 CASE( 2 ) ! scaling with mean sea-ice thickness 467 DO_2D_00_00 468 #if defined key_si3 469 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * hm_i(ji,jj) * 2. ) * tmask(ji,jj,1) 470 #elif defined key_cice 471 zmaxice = MAXVAL( h_i(ji,jj,:) ) 472 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 473 #endif 474 END_2D 475 ! 476 CASE( 3 ) ! scaling with max sea-ice thickness 477 DO_2D_00_00 478 zmaxice = MAXVAL( h_i(ji,jj,:) ) 479 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 480 END_2D 481 ! 482 END SELECT 483 #endif 484 ! 485 DO_2D_00_00 486 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 487 END_2D 488 ! 489 ELSE 447 490 zmxlm(:,:,1) = rn_mxl0 448 491 ENDIF 492 449 493 ! 450 494 DO_3D_00_00( 2, jpkm1 ) … … 518 562 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 519 563 DO_3D_00_00( 2, jpkm1 ) 520 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)564 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 521 565 END_3D 522 566 ENDIF … … 550 594 INTEGER :: ios 551 595 !! 552 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , & 553 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 554 & rn_mxl0 , nn_pdl , ln_drg , ln_lc , rn_lc, & 555 & nn_etau , nn_htau , rn_efr , rn_eice 596 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , & 597 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 598 & rn_mxl0 , nn_mxlice, rn_mxlice, & 599 & nn_pdl , ln_drg , ln_lc , rn_lc, & 600 & nn_etau , nn_htau , rn_efr , rn_eice 556 601 !!---------------------------------------------------------------------- 557 602 ! … … 579 624 WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl 580 625 WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0 626 IF( ln_mxl0 ) THEN 627 WRITE(numout,*) ' type of scaling under sea-ice nn_mxlice = ', nn_mxlice 628 IF( nn_mxlice == 1 ) & 629 WRITE(numout,*) ' ice thickness when scaling under sea-ice rn_mxlice = ', rn_mxlice 630 ENDIF 581 631 WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0 582 632 WRITE(numout,*) ' top/bottom friction forcing flag ln_drg = ', ln_drg -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/nemogcm.F90
r12489 r13159 84 84 #endif 85 85 ! 86 USE in_out_manager ! I/O manager 86 87 USE lib_mpp ! distributed memory computing 87 88 USE mppini ! shared/distributed memory setting (mpp_init routine) … … 185 186 END DO 186 187 ! 187 IF( .NOT. Agrif_Root() ) THEN188 CALL Agrif_ParentGrid_To_ChildGrid()189 IF( ln_diaobs ) CALL dia_obs_wri190 IF( ln_timing ) CALL timing_finalize191 CALL Agrif_ChildGrid_To_ParentGrid()192 ENDIF193 !194 188 # else 195 189 ! … … 236 230 IF( nstop /= 0 .AND. lwp ) THEN ! error print 237 231 WRITE(ctmp1,*) ' ==>>> nemo_gcm: a total of ', nstop, ' errors have been found' 238 CALL ctl_stop( ctmp1 ) 232 IF( ngrdstop > 0 ) THEN 233 WRITE(ctmp9,'(i2)') ngrdstop 234 WRITE(ctmp2,*) ' E R R O R detected in Agrif grid '//TRIM(ctmp9) 235 WRITE(ctmp3,*) ' Look for "E R R O R" messages in all existing '//TRIM(ctmp9)//'_ocean_output* files' 236 CALL ctl_stop( ' ', ctmp1, ' ', ctmp2, ' ', ctmp3 ) 237 ELSE 238 WRITE(ctmp2,*) ' Look for "E R R O R" messages in all existing ocean_output* files' 239 CALL ctl_stop( ' ', ctmp1, ' ', ctmp2 ) 240 ENDIF 239 241 ENDIF 240 242 ! … … 248 250 #else 249 251 IF ( lk_oasis ) THEN ; CALL cpl_finalize ! end coupling and mpp communications with OASIS 250 ELSEIF( lk_mpp ) THEN ; CALL mppstop ! end mpp communications252 ELSEIF( lk_mpp ) THEN ; CALL mppstop ! end mpp communications 251 253 ENDIF 252 254 #endif … … 317 319 IF( lwm ) CALL ctl_opn( numond, 'output.namelist.dyn', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. ) 318 320 ! open /dev/null file to be able to supress output write easily 321 IF( Agrif_Root() ) THEN 319 322 CALL ctl_opn( numnul, '/dev/null', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. ) 320 ! 323 #ifdef key_agrif 324 ELSE 325 numnul = Agrif_Parent(numnul) 326 #endif 327 ENDIF 321 328 ! !--------------------! 322 329 ! ! Open listing units ! -> need sn_cfctl from namctl to define lwp … … 329 336 ! 330 337 ! finalize the definition of namctl variables 331 IF( sn_cfctl%l_allon ) THEN 332 ! Turn on all options. 333 CALL nemo_set_cfctl( sn_cfctl, .TRUE., .TRUE. ) 334 ! Ensure all processors are active 335 sn_cfctl%procmin = 0 ; sn_cfctl%procmax = 1000000 ; sn_cfctl%procincr = 1 336 ELSEIF( sn_cfctl%l_config ) THEN 337 ! Activate finer control of report outputs 338 ! optionally switch off output from selected areas (note this only 339 ! applies to output which does not involve global communications) 340 IF( ( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax ) .OR. & 341 & ( MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 ) ) & 342 & CALL nemo_set_cfctl( sn_cfctl, .FALSE., .FALSE. ) 343 ELSE 344 ! turn off all options. 345 CALL nemo_set_cfctl( sn_cfctl, .FALSE., .TRUE. ) 346 ENDIF 338 IF( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax .OR. MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 ) & 339 & CALL nemo_set_cfctl( sn_cfctl, .FALSE. ) 347 340 ! 348 341 lwp = (narea == 1) .OR. sn_cfctl%l_oceout ! control of all listing output print … … 528 521 WRITE(numout,*) '~~~~~~~~' 529 522 WRITE(numout,*) ' Namelist namctl' 530 WRITE(numout,*) ' sn_cfctl%l_glochk = ', sn_cfctl%l_glochk531 WRITE(numout,*) ' sn_cfctl%l_allon = ', sn_cfctl%l_allon532 WRITE(numout,*) ' finer control over o/p sn_cfctl%l_config = ', sn_cfctl%l_config533 523 WRITE(numout,*) ' sn_cfctl%l_runstat = ', sn_cfctl%l_runstat 534 524 WRITE(numout,*) ' sn_cfctl%l_trcstat = ', sn_cfctl%l_trcstat … … 678 668 679 669 680 SUBROUTINE nemo_set_cfctl(sn_cfctl, setto , for_all)670 SUBROUTINE nemo_set_cfctl(sn_cfctl, setto ) 681 671 !!---------------------------------------------------------------------- 682 672 !! *** ROUTINE nemo_set_cfctl *** 683 673 !! 684 674 !! ** Purpose : Set elements of the output control structure to setto. 685 !! for_all should be .false. unless all areas are to be686 !! treated identically.687 675 !! 688 676 !! ** Method : Note this routine can be used to switch on/off some 689 !! types of output for selected areas but any output types 690 !! that involve global communications (e.g. mpp_max, glob_sum) 691 !! should be protected from selective switching by the 692 !! for_all argument 693 !!---------------------------------------------------------------------- 694 LOGICAL :: setto, for_all 695 TYPE(sn_ctl) :: sn_cfctl 696 !!---------------------------------------------------------------------- 697 IF( for_all ) THEN 698 sn_cfctl%l_runstat = setto 699 sn_cfctl%l_trcstat = setto 700 ENDIF 677 !! types of output for selected areas. 678 !!---------------------------------------------------------------------- 679 TYPE(sn_ctl), INTENT(inout) :: sn_cfctl 680 LOGICAL , INTENT(in ) :: setto 681 !!---------------------------------------------------------------------- 682 sn_cfctl%l_runstat = setto 683 sn_cfctl%l_trcstat = setto 701 684 sn_cfctl%l_oceout = setto 702 685 sn_cfctl%l_layout = setto -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/step.F90
r12489 r13159 82 82 !!---------------------------------------------------------------------- 83 83 INTEGER :: ji, jj, jk ! dummy loop indice 84 INTEGER :: indic ! error indicator if < 085 84 !!gm kcall can be removed, I guess 86 85 INTEGER :: kcall ! optional integer argument (dom_vvl_sf_nxt) 87 86 !! --------------------------------------------------------------------- 88 87 #if defined key_agrif 88 IF( nstop > 0 ) RETURN ! avoid to go further if an error was detected during previous time step (child grid) 89 89 kstp = nit000 + Agrif_Nb_Step() 90 90 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices … … 114 114 ! update I/O and calendar 115 115 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 116 indic = 0 ! reset to no error condition117 118 116 IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS) 119 CALL iom_init( cxios_context, ld_closedef=.FALSE. ) ! for model grid (including p assible AGRIF zoom)117 CALL iom_init( cxios_context, ld_closedef=.FALSE. ) ! for model grid (including possible AGRIF zoom) 120 118 IF( lk_diamlr ) CALL dia_mlr_iom_init ! with additional setup for multiple-linear-regression analysis 121 119 CALL iom_init_closedef … … 309 307 #if defined key_agrif 310 308 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 311 ! AGRIF 309 ! AGRIF recursive integration 312 310 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 313 311 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices 314 312 CALL Agrif_Integrate_ChildGrids( stp ) ! allows to finish all the Child Grids before updating 315 313 316 IF( Agrif_NbStepint() == 0 ) THEN 317 CALL Agrif_update_all( ) ! Update all components 318 ENDIF 319 #endif 320 IF( ln_diaobs ) CALL dia_obs ( kstp, Nnn ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 321 314 #endif 322 315 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 323 316 ! Control 324 317 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 325 CALL stp_ctl ( kstp, Nbb, Nnn, indic ) 326 318 CALL stp_ctl ( kstp, Nnn ) 319 320 #if defined key_agrif 321 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 322 ! AGRIF update 323 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 324 IF( Agrif_NbStepint() == 0 .AND. nstop == 0 ) THEN 325 CALL Agrif_update_all( ) ! Update all components 326 ENDIF 327 328 #endif 329 IF( ln_diaobs .AND. nstop == 0 ) CALL dia_obs( kstp, Nnn ) ! obs-minus-model (assimilation) diags (after dynamics update) 330 331 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 332 ! File manipulation at the end of the first time step 333 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 327 334 IF( kstp == nit000 ) THEN ! 1st time step only 328 335 CALL iom_close( numror ) ! close input ocean restart file … … 334 341 ! Coupled mode 335 342 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 336 !!gm why lk_oasis and not lk_cpl ???? 337 IF( lk_oasis ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges 343 IF( lk_oasis .AND. nstop == 0 ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges 338 344 ! 339 345 #if defined key_iomput 340 IF( kstp == nitend .OR. indic < 0 ) THEN 346 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 347 ! Finalize contextes if end of simulation or error detected 348 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 349 IF( kstp == nitend .OR. nstop > 0 ) THEN 341 350 CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF 342 IF(lrxios) CALL iom_context_finalize( crxios_context)351 IF( lrxios ) CALL iom_context_finalize( crxios_context ) 343 352 IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 344 353 ENDIF -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/stpctl.F90
r12377 r13159 19 19 USE dom_oce ! ocean space and time domain variables 20 20 USE c1d ! 1D vertical configuration 21 USE zdf_oce , ONLY : ln_zad_Aimp ! ocean vertical physics variables 22 USE wet_dry, ONLY : ll_wd, ssh_ref ! reference depth for negative bathy 23 ! 21 24 USE diawri ! Standard run outputs (dia_wri_state routine) 22 !23 25 USE in_out_manager ! I/O manager 24 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 25 27 USE lib_mpp ! distributed memory computing 26 USE zdf_oce , ONLY : ln_zad_Aimp ! ocean vertical physics variables 27 USE wet_dry, ONLY : ll_wd, ssh_ref ! reference depth for negative bathy 28 28 ! 29 29 USE netcdf ! NetCDF library 30 30 IMPLICIT NONE … … 33 33 PUBLIC stp_ctl ! routine called by step.F90 34 34 35 INTEGER :: idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus36 LOGICAL :: lsomeoce35 INTEGER :: nrunid ! netcdf file id 36 INTEGER, DIMENSION(8) :: nvarid ! netcdf variable id 37 37 !!---------------------------------------------------------------------- 38 38 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 42 42 CONTAINS 43 43 44 SUBROUTINE stp_ctl( kt, K bb, Kmm, kindic)44 SUBROUTINE stp_ctl( kt, Kmm ) 45 45 !!---------------------------------------------------------------------- 46 46 !! *** ROUTINE stp_ctl *** … … 50 50 !! ** Method : - Save the time step in numstp 51 51 !! - Print it each 50 time steps 52 !! - Stop the run IF problem encountered by setting indic=-352 !! - Stop the run IF problem encountered by setting nstop > 0 53 53 !! Problems checked: |ssh| maximum larger than 10 m 54 54 !! |U| maximum larger than 10 m/s … … 57 57 !! ** Actions : "time.step" file = last ocean time-step 58 58 !! "run.stat" file = run statistics 59 !! nstop indicator sheared among all local domain (lk_mpp=T)59 !! nstop indicator sheared among all local domain 60 60 !!---------------------------------------------------------------------- 61 61 INTEGER, INTENT(in ) :: kt ! ocean time-step index 62 INTEGER, INTENT(in ) :: Kbb, Kmm ! ocean time level index 63 INTEGER, INTENT(inout) :: kindic ! error indicator 64 !! 65 INTEGER :: ji, jj, jk ! dummy loop indices 66 INTEGER, DIMENSION(2) :: ih ! min/max loc indices 67 INTEGER, DIMENSION(3) :: iu, is1, is2 ! min/max loc indices 68 REAL(wp) :: zzz ! local real 69 REAL(wp), DIMENSION(9) :: zmax 70 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 71 CHARACTER(len=20) :: clname 72 !!---------------------------------------------------------------------- 73 ! 74 ll_wrtstp = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 75 ll_colruns = ll_wrtstp .AND. ( sn_cfctl%l_runstat ) 76 ll_wrtruns = ll_colruns .AND. lwm 77 IF( kt == nit000 .AND. lwp ) THEN 78 WRITE(numout,*) 79 WRITE(numout,*) 'stp_ctl : time-stepping control' 80 WRITE(numout,*) '~~~~~~~' 81 ! ! open time.step file 82 IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 83 ! ! open run.stat file(s) at start whatever 84 ! ! the value of sn_cfctl%ptimincr 85 IF( lwm .AND. ( sn_cfctl%l_runstat ) ) THEN 62 INTEGER, INTENT(in ) :: Kmm ! ocean time level index 63 !! 64 INTEGER :: ji ! dummy loop indices 65 INTEGER :: idtime, istatus 66 INTEGER , DIMENSION(9) :: iareasum, iareamin, iareamax 67 INTEGER , DIMENSION(3,4) :: iloc ! min/max loc indices 68 REAL(wp) :: zzz ! local real 69 REAL(wp), DIMENSION(9) :: zmax, zmaxlocal 70 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 71 LOGICAL, DIMENSION(jpi,jpj,jpk) :: llmsk 72 CHARACTER(len=20) :: clname 73 !!---------------------------------------------------------------------- 74 IF( nstop > 0 .AND. ngrdstop > -1 ) RETURN ! stpctl was already called by a child grid 75 ! 76 ll_wrtstp = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 77 ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1 78 ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm 79 ! 80 IF( kt == nit000 ) THEN 81 ! 82 IF( lwp ) THEN 83 WRITE(numout,*) 84 WRITE(numout,*) 'stp_ctl : time-stepping control' 85 WRITE(numout,*) '~~~~~~~' 86 ENDIF 87 ! ! open time.step ascii file, done only by 1st subdomain 88 IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 89 ! 90 IF( ll_wrtruns ) THEN 91 ! ! open run.stat ascii file, done only by 1st subdomain 86 92 CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 93 ! ! open run.stat.nc netcdf file, done only by 1st subdomain 87 94 clname = 'run.stat.nc' 88 95 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 89 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, idrun)90 istatus = NF90_DEF_DIM( idrun, 'time', NF90_UNLIMITED, idtime )91 istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh)92 istatus = NF90_DEF_VAR( idrun, 'abs_u_max', NF90_DOUBLE, (/ idtime /), idu)93 istatus = NF90_DEF_VAR( idrun, 's_min', NF90_DOUBLE, (/ idtime /), ids1)94 istatus = NF90_DEF_VAR( idrun, 's_max', NF90_DOUBLE, (/ idtime /), ids2)95 istatus = NF90_DEF_VAR( idrun, 't_min', NF90_DOUBLE, (/ idtime /), idt1)96 istatus = NF90_DEF_VAR( idrun, 't_max', NF90_DOUBLE, (/ idtime /), idt2)96 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid ) 97 istatus = NF90_DEF_DIM( nrunid, 'time', NF90_UNLIMITED, idtime ) 98 istatus = NF90_DEF_VAR( nrunid, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), nvarid(1) ) 99 istatus = NF90_DEF_VAR( nrunid, 'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2) ) 100 istatus = NF90_DEF_VAR( nrunid, 's_min', NF90_DOUBLE, (/ idtime /), nvarid(3) ) 101 istatus = NF90_DEF_VAR( nrunid, 's_max', NF90_DOUBLE, (/ idtime /), nvarid(4) ) 102 istatus = NF90_DEF_VAR( nrunid, 't_min', NF90_DOUBLE, (/ idtime /), nvarid(5) ) 103 istatus = NF90_DEF_VAR( nrunid, 't_max', NF90_DOUBLE, (/ idtime /), nvarid(6) ) 97 104 IF( ln_zad_Aimp ) THEN 98 istatus = NF90_DEF_VAR( idrun, 'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1)99 istatus = NF90_DEF_VAR( idrun, 'Cf_max', NF90_DOUBLE, (/ idtime /), idc1)105 istatus = NF90_DEF_VAR( nrunid, 'Cf_max', NF90_DOUBLE, (/ idtime /), nvarid(7) ) 106 istatus = NF90_DEF_VAR( nrunid,'abs_wi_max',NF90_DOUBLE, (/ idtime /), nvarid(8) ) 100 107 ENDIF 101 istatus = NF90_ENDDEF(idrun) 102 zmax(8:9) = 0._wp ! initialise to zero in case ln_zad_Aimp option is not in use 103 ENDIF 104 ENDIF 105 IF( kt == nit000 ) lsomeoce = COUNT( ssmask(:,:) == 1._wp ) > 0 106 ! 107 IF(lwm .AND. ll_wrtstp) THEN !== current time step ==! ("time.step" file) 108 istatus = NF90_ENDDEF(nrunid) 109 ENDIF 110 ! 111 ENDIF 112 ! 113 ! !== write current time step ==! 114 ! !== done only by 1st subdomain at writting timestep ==! 115 IF( lwm .AND. ll_wrtstp ) THEN 108 116 WRITE ( numstp, '(1x, i8)' ) kt 109 117 REWIND( numstp ) 110 118 ENDIF 111 ! 112 ! !== test of extrema ==! 113 IF( ll_wd ) THEN 114 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) ) ) ! ssh max 115 ELSE 116 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) ) ) ! ssh max 117 ENDIF 118 zmax(2) = MAXVAL( ABS( uu(:,:,:,Kmm) ) ) ! velocity max (zonal only) 119 zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! minus salinity max 120 zmax(4) = MAXVAL( ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! salinity max 121 zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! minus temperature max 122 zmax(6) = MAXVAL( ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! temperature max 123 zmax(7) = REAL( nstop , wp ) ! stop indicator 124 IF( ln_zad_Aimp ) THEN 125 zmax(8) = MAXVAL( ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 126 zmax(9) = MAXVAL( Cu_adv(:,:,:) , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max 127 ENDIF 128 ! 119 ! !== test of local extrema ==! 120 ! !== done by all processes at every time step ==! 121 ! 122 ! define zmax default value. needed for land processors 123 IF( ll_colruns ) THEN ! default value: must not be kept when calling mpp_max -> must be as small as possible 124 zmax(:) = -HUGE(1._wp) 125 ELSE ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...) 126 zmax(:) = 0._wp 127 zmax(3) = -1._wp ! avoid salinity minimum at 0. 128 ENDIF 129 ! 130 llmsk(:,:,1) = ssmask(:,:) == 1._wp 131 IF( COUNT( llmsk(:,:,1) ) > 0 ) THEN ! avoid huge values sent back for land processors... 132 IF( ll_wd ) THEN 133 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) ) ! ssh max 134 ELSE 135 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) ), mask = llmsk(:,:,1) ) ! ssh max 136 ENDIF 137 ENDIF 138 zmax(2) = MAXVAL( ABS( uu(:,:,:,Kmm) ) ) ! velocity max (zonal only) 139 llmsk(:,:,:) = tmask(:,:,:) == 1._wp 140 IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN ! avoid huge values sent back for land processors... 141 zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk ) ! minus salinity max 142 zmax(4) = MAXVAL( ts(:,:,:,jp_sal,Kmm), mask = llmsk ) ! salinity max 143 IF( ll_colruns .OR. jpnij == 1 ) THEN ! following variables are used only in the netcdf file 144 zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk ) ! minus temperature max 145 zmax(6) = MAXVAL( ts(:,:,:,jp_tem,Kmm), mask = llmsk ) ! temperature max 146 IF( ln_zad_Aimp ) THEN 147 zmax(7) = MAXVAL( Cu_adv(:,:,:) , mask = llmsk ) ! partitioning coeff. max 148 llmsk(:,:,:) = wmask(:,:,:) == 1._wp 149 IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN ! avoid huge values sent back for land processors... 150 zmax(8) = MAXVAL(ABS( wi(:,:,:) ), mask = llmsk ) ! implicit vertical vel. max 151 ENDIF 152 ENDIF 153 ENDIF 154 ENDIF 155 zmax(9) = REAL( nstop, wp ) ! stop indicator 156 ! !== get global extrema ==! 157 ! !== done by all processes if writting run.stat ==! 129 158 IF( ll_colruns ) THEN 159 zmaxlocal(:) = zmax(:) 130 160 CALL mpp_max( "stpctl", zmax ) ! max over the global domain 131 nstop = NINT( zmax(7) ) ! nstop indicator sheared among all local domains 132 ENDIF 133 ! !== run statistics ==! ("run.stat" files) 161 nstop = NINT( zmax(9) ) ! update nstop indicator (now sheared among all local domains) 162 ENDIF 163 ! !== write "run.stat" files ==! 164 ! !== done only by 1st subdomain at writting timestep ==! 134 165 IF( ll_wrtruns ) THEN 135 166 WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 136 istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) )137 istatus = NF90_PUT_VAR( idrun, idu, (/ zmax(2)/), (/kt/), (/1/) )138 istatus = NF90_PUT_VAR( idrun, ids1, (/-zmax(3)/), (/kt/), (/1/) )139 istatus = NF90_PUT_VAR( idrun, ids2, (/ zmax(4)/), (/kt/), (/1/) )140 istatus = NF90_PUT_VAR( idrun, idt1, (/-zmax(5)/), (/kt/), (/1/) )141 istatus = NF90_PUT_VAR( idrun, idt2, (/ zmax(6)/), (/kt/), (/1/) )167 istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 168 istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 169 istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) ) 170 istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) ) 171 istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) ) 172 istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) ) 142 173 IF( ln_zad_Aimp ) THEN 143 istatus = NF90_PUT_VAR( idrun, idw1, (/ zmax(8)/), (/kt/), (/1/) ) 144 istatus = NF90_PUT_VAR( idrun, idc1, (/ zmax(9)/), (/kt/), (/1/) ) 145 ENDIF 146 IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun) 147 IF( kt == nitend ) istatus = NF90_CLOSE(idrun) 174 istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) ) 175 istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) ) 176 ENDIF 177 IF( kt == nitend ) istatus = NF90_CLOSE(nrunid) 148 178 END IF 149 ! !== error handling ==! 150 IF( ( sn_cfctl%l_glochk .OR. lsomeoce ) .AND. ( & ! domain contains some ocean points, check for sensible ranges 151 & zmax(1) > 20._wp .OR. & ! too large sea surface height ( > 20 m ) 152 & zmax(2) > 10._wp .OR. & ! too large velocity ( > 10 m/s) 153 & zmax(3) >= 0._wp .OR. & ! negative or zero sea surface salinity 154 & zmax(4) >= 100._wp .OR. & ! too large sea surface salinity ( > 100 ) 155 & zmax(4) < 0._wp .OR. & ! too large sea surface salinity (keep this line for sea-ice) 156 & ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN ! NaN encounter in the tests 157 IF( lk_mpp .AND. sn_cfctl%l_glochk ) THEN 158 ! have use mpp_max (because sn_cfctl%l_glochk=.T. and distributed) 159 CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm)) , ssmask(:,:) , zzz, ih ) 160 CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm)) , umask (:,:,:), zzz, iu ) 161 CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 ) 162 CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 ) 179 ! !== error handling ==! 180 ! !== done by all processes at every time step ==! 181 ! 182 IF( zmax(1) > 20._wp .OR. & ! too large sea surface height ( > 20 m ) 183 & zmax(2) > 10._wp .OR. & ! too large velocity ( > 10 m/s) 184 & zmax(3) >= 0._wp .OR. & ! negative or zero sea surface salinity 185 & zmax(4) >= 100._wp .OR. & ! too large sea surface salinity ( > 100 ) 186 & zmax(4) < 0._wp .OR. & ! too large sea surface salinity (keep this line for sea-ice) 187 & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR. & ! NaN encounter in the tests 188 & ABS( zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN ! Infinity encounter in the tests 189 ! 190 iloc(:,:) = 0 191 IF( ll_colruns ) THEN ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc 192 ! first: close the netcdf file, so we can read it 193 IF( lwm .AND. kt /= nitend ) istatus = NF90_CLOSE(nrunid) 194 ! get global loc on the min/max 195 CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:, Kmm)), ssmask(:,: ), zzz, iloc(1:2,1) ) ! mpp_maxloc ok if mask = F 196 CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:, Kmm)), umask(:,:,:), zzz, iloc(1:3,2) ) 197 CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm) , tmask(:,:,:), zzz, iloc(1:3,3) ) 198 CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm) , tmask(:,:,:), zzz, iloc(1:3,4) ) 199 ! find which subdomain has the max. 200 iareamin(:) = jpnij+1 ; iareamax(:) = 0 ; iareasum(:) = 0 201 DO ji = 1, 9 202 IF( zmaxlocal(ji) == zmax(ji) ) THEN 203 iareamin(ji) = narea ; iareamax(ji) = narea ; iareasum(ji) = 1 204 ENDIF 205 END DO 206 CALL mpp_min( "stpctl", iareamin ) ! min over the global domain 207 CALL mpp_max( "stpctl", iareamax ) ! max over the global domain 208 CALL mpp_sum( "stpctl", iareasum ) ! sum over the global domain 209 ELSE ! find local min and max locations: 210 ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc 211 iloc(1:2,1) = MAXLOC( ABS( ssh(:,:, Kmm)), mask = ssmask(:,: ) == 1._wp ) + (/ nimpp - 1, njmpp - 1 /) 212 iloc(1:3,2) = MAXLOC( ABS( uu(:,:,:, Kmm)), mask = umask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 213 iloc(1:3,3) = MINLOC( ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 214 iloc(1:3,4) = MAXLOC( ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 215 iareamin(:) = narea ; iareamax(:) = narea ; iareasum(:) = 1 ! this is local information 216 ENDIF 217 ! 218 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests' 219 CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 220 CALL wrt_line( ctmp3, kt, '|U| max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 221 CALL wrt_line( ctmp4, kt, 'Sal min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 222 CALL wrt_line( ctmp5, kt, 'Sal max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 223 IF( Agrif_Root() ) THEN 224 WRITE(ctmp6,*) ' ===> output of last computed fields in output.abort* files' 163 225 ELSE 164 ! find local min and max locations 165 ih(:) = MAXLOC( ABS( ssh(:,:,Kmm) ) ) + (/ nimpp - 1, njmpp - 1 /) 166 iu(:) = MAXLOC( ABS( uu (:,:,:,Kmm) ) ) + (/ nimpp - 1, njmpp - 1, 0 /) 167 is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 168 is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 169 ENDIF 170 171 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests' 172 WRITE(ctmp2,9100) kt, zmax(1), ih(1) , ih(2) 173 WRITE(ctmp3,9200) kt, zmax(2), iu(1) , iu(2) , iu(3) 174 WRITE(ctmp4,9300) kt, - zmax(3), is1(1), is1(2), is1(3) 175 WRITE(ctmp5,9400) kt, zmax(4), is2(1), is2(2), is2(3) 176 WRITE(ctmp6,*) ' ===> output of last computed fields in output.abort.nc file' 177 226 WRITE(ctmp6,*) ' ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files' 227 ENDIF 228 ! 178 229 CALL dia_wri_state( Kmm, 'output.abort' ) ! create an output.abort file 179 180 IF( .NOT. sn_cfctl%l_glochk ) THEN181 WRITE(ctmp8,*) 'E R R O R message from sub-domain: ', narea182 CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ctmp6)183 ELSE184 CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6, ' ' )185 ENDIF186 187 kindic = -3188 !189 ENDIF190 !191 9100 FORMAT (' kt=',i8,' |ssh| max: ',1pg11.4,', at i j : ',2i5) 192 9200 FORMAT (' kt=',i8,' |U| max: ',1pg11.4,', at i j k: ',3i5) 193 9300 FORMAT (' kt=',i8,' S min: ',1pg11.4,', at i j k: ',3i5) 194 9400 FORMAT (' kt=',i8,' S max: ',1pg11.4,', at i j k: ',3i5) 230 ! 231 IF( ll_colruns .or. jpnij == 1 ) THEN ! all processes synchronized -> use lwp to print in opened ocean.output files 232 IF(lwp) THEN ; CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 ) 233 ELSE ; nstop = MAX(1, nstop) ! make sure nstop > 0 (automatically done when calling ctl_stop) 234 ENDIF 235 ELSE ! only mpi subdomains with errors are here -> STOP now 236 CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 ) 237 ENDIF 238 ! 239 ENDIF 240 ! 241 IF( nstop > 0 ) THEN ! an error was detected and we did not abort yet... 242 ngrdstop = Agrif_Fixed() ! store which grid got this error 243 IF( .NOT. ll_colruns .AND. jpnij > 1 ) CALL ctl_stop( 'STOP' ) ! we must abort here to avoid MPI deadlock 244 ENDIF 245 ! 195 246 9500 FORMAT(' it :', i8, ' |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 196 247 ! 197 248 END SUBROUTINE stp_ctl 249 250 251 SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax ) 252 !!---------------------------------------------------------------------- 253 !! *** ROUTINE wrt_line *** 254 !! 255 !! ** Purpose : write information line 256 !! 257 !!---------------------------------------------------------------------- 258 CHARACTER(len=*), INTENT( out) :: cdline 259 CHARACTER(len=*), INTENT(in ) :: cdprefix 260 REAL(wp), INTENT(in ) :: pval 261 INTEGER, DIMENSION(3), INTENT(in ) :: kloc 262 INTEGER, INTENT(in ) :: kt, ksum, kmin, kmax 263 ! 264 CHARACTER(len=80) :: clsuff 265 CHARACTER(len=9 ) :: clkt, clsum, clmin, clmax 266 CHARACTER(len=9 ) :: cli, clj, clk 267 CHARACTER(len=1 ) :: clfmt 268 CHARACTER(len=4 ) :: cl4 ! needed to be able to compile with Agrif, I don't know why 269 INTEGER :: ifmtk 270 !!---------------------------------------------------------------------- 271 WRITE(clkt , '(i9)') kt 272 273 WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij ,wp))) + 1 ! how many digits to we need to write ? (we decide max = 9) 274 !!! WRITE(clsum, '(i'//clfmt//')') ksum ! this is creating a compilation error with AGRIF 275 cl4 = '(i'//clfmt//')' ; WRITE(clsum, cl4) ksum 276 WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1 ! how many digits to we need to write ? (we decide max = 9) 277 cl4 = '(i'//clfmt//')' ; WRITE(clmin, cl4) kmin-1 278 WRITE(clmax, cl4) kmax-1 279 ! 280 WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1 ! how many digits to we need to write jpiglo? (we decide max = 9) 281 cl4 = '(i'//clfmt//')' ; WRITE(cli, cl4) kloc(1) ! this is ok with AGRIF 282 WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1 ! how many digits to we need to write jpjglo? (we decide max = 9) 283 cl4 = '(i'//clfmt//')' ; WRITE(clj, cl4) kloc(2) ! this is ok with AGRIF 284 ! 285 IF( ksum == 1 ) THEN ; WRITE(clsuff,9100) TRIM(clmin) 286 ELSE ; WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax) 287 ENDIF 288 IF(kloc(3) == 0) THEN 289 ifmtk = INT(LOG10(REAL(jpk,wp))) + 1 ! how many digits to we need to write jpk? (we decide max = 9) 290 clk = REPEAT(' ', ifmtk) ! create the equivalent in blank string 291 WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff) 292 ELSE 293 WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1 ! how many digits to we need to write jpk? (we decide max = 9) 294 !!! WRITE(clk, '(i'//clfmt//')') kloc(3) ! this is creating a compilation error with AGRIF 295 cl4 = '(i'//clfmt//')' ; WRITE(clk, cl4) kloc(3) ! this is ok with AGRIF 296 WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), TRIM(clk), TRIM(clsuff) 297 ENDIF 298 ! 299 9100 FORMAT('MPI rank ', a) 300 9200 FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a) 301 9300 FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j ', a, ' ', a, ' ', a, ' ', a) 302 9400 FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a) 303 ! 304 END SUBROUTINE wrt_line 305 198 306 199 307 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.