- Timestamp:
- 2015-07-20T19:43:15+02:00 (9 years ago)
- Location:
- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r5123 r5619 43 43 INTEGER , PUBLIC :: nn_closea !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 44 44 INTEGER , PUBLIC :: nn_euler !: =0 start with forward time step or not (=1) 45 LOGICAL , PUBLIC :: ln_iscpl !: coupling with ice sheet 45 46 LOGICAL , PUBLIC :: ln_crs !: Apply grid coarsening to dynamical model output or online passive tracers 46 47 … … 252 253 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level 253 254 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku, mbkv !: vertical index of the bottom last U- and W- ocean level 254 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy 255 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i , umask_i, vmask_i, fmask_i!: interior domain T-point mask256 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask 255 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 256 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask 257 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function 257 258 258 259 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfdep !: top first ocean level (ISF) 259 260 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: first wet T-, U-, V-, F- ocean level (ISF) 260 261 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfdep !: Iceshelf draft (ISF) 261 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask !: surface domain T-point mask 262 262 263 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts 263 264 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts 264 265 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask !: land/ocean mask at WT-, WU- and WV-pts … … 389 390 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 390 391 391 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & 392 & tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 393 & bmask(jpi,jpj) , & 394 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 392 ALLOCATE( mbathy(jpi,jpj) , bathy (jpi,jpj) , & 393 & tmask_i(jpi,jpj) , & 394 & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 395 & bmask(jpi,jpj) , & 396 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(9) ) 395 397 396 398 ! (ISF) Allocation of basic array 397 399 ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), & 398 400 & mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) , & 399 & mikf(jpi,jpj), ssmask(jpi,jpj),STAT=ierr(10) )401 & mikf(jpi,jpj), STAT=ierr(10) ) 400 402 401 403 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk), & -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r5363 r5619 111 111 END DO 112 112 ! ! Inverse of the local depth 113 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)114 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)113 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 114 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 115 115 116 116 CALL dom_stp ! time step 117 IF( nmsh /= 0 117 IF( nmsh /= 0 .AND. .NOT. ln_iscpl ) CALL dom_wri ! Create a domain file 118 118 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 119 119 ! … … 138 138 & nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl, & 139 139 & nn_it000, nn_itend , nn_date0 , nn_leapy , nn_istate , nn_stock , & 140 & nn_write, ln_dimgnnn, ln_mskland , ln_cfmeta , ln_clobber, nn_chunksz, nn_euler 140 & nn_write, ln_dimgnnn, ln_mskland , ln_cfmeta , ln_clobber, nn_chunksz, nn_euler, ln_iscpl 141 141 NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, & 142 142 & nn_acc , rn_atfp , rn_rdt , rn_rdtmin , & … … 192 192 WRITE(numout,*) ' overwrite an existing file ln_clobber = ', ln_clobber 193 193 WRITE(numout,*) ' NetCDF chunksize (bytes) nn_chunksz = ', nn_chunksz 194 WRITE(numout,*) ' IS coupling at the restart step ln_iscpl = ', ln_iscpl 194 195 ENDIF 195 196 -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r5552 r5619 264 264 END DO 265 265 END DO 266 ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point266 ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point 267 267 DO jj = 1, jpjm1 268 268 DO ji = 1, fs_jpim1 ! vector loop 269 umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:)))270 vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:)))269 ssumask(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:))) 270 ssvmask(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:))) 271 271 END DO 272 272 DO ji = 1, jpim1 ! NO vector opt. 273 fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) &273 ssfmask(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) & 274 274 & * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 275 275 END DO 276 276 END DO 277 CALL lbc_lnk( umask , 'U', 1._wp ) ! Lateral boundary conditions278 CALL lbc_lnk( vmask , 'V', 1._wp )279 CALL lbc_lnk( fmask , 'F', 1._wp )280 CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions281 CALL lbc_lnk( vmask_i, 'V', 1._wp )282 CALL lbc_lnk( fmask_i, 'F', 1._wp )277 CALL lbc_lnk( umask , 'U', 1._wp ) ! Lateral boundary conditions 278 CALL lbc_lnk( vmask , 'V', 1._wp ) 279 CALL lbc_lnk( fmask , 'F', 1._wp ) 280 CALL lbc_lnk( ssumask, 'U', 1._wp ) ! Lateral boundary conditions 281 CALL lbc_lnk( ssvmask, 'V', 1._wp ) 282 CALL lbc_lnk( ssfmask, 'F', 1._wp ) 283 283 284 284 ! 3. Ocean/land mask at wu-, wv- and w points -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
r3294 r5619 28 28 CONTAINS 29 29 30 SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid )30 SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, kkk ) 31 31 !!---------------------------------------------------------------------- 32 32 !! *** ROUTINE dom_ngb *** … … 40 40 REAL(wp) , INTENT(in ) :: plon, plat ! longitude,latitude of the point 41 41 INTEGER , INTENT( out) :: kii, kjj ! i-,j-index of the closes grid point 42 INTEGER , INTENT(in ), OPTIONAL :: kkk 42 43 CHARACTER(len=1), INTENT(in ) :: cdgrid ! grid name 'T', 'U', 'V', 'W' 43 44 ! 44 45 INTEGER , DIMENSION(2) :: iloc 46 INTEGER :: jk 45 47 REAL(wp) :: zlon, zmini 46 48 REAL(wp), POINTER, DIMENSION(:,:) :: zglam, zgphi, zmask, zdist … … 52 54 ! 53 55 zmask(:,:) = 0._wp 56 jk = 1 57 IF (PRESENT(kkk)) jk=kkk 54 58 SELECT CASE( cdgrid ) 55 CASE( 'U' ) ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej, 1)56 CASE( 'V' ) ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej, 1)57 CASE( 'F' ) ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej, 1)58 CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej, 1)59 CASE( 'U' ) ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,jk) 60 CASE( 'V' ) ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,jk) 61 CASE( 'F' ) ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,jk) 62 CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,jk) 59 63 END SELECT 60 64 61 zlon = MOD( plon + 720., 360. ) ! plon between 0 and 360 62 zglam(:,:) = MOD( zglam(:,:) + 720., 360. ) ! glam between 0 and 360 63 IF( zlon > 270. ) zlon = zlon - 360. ! zlon between -90 and 270 64 IF( zlon < 90. ) WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360. ! glam between -180 and 180 65 66 zglam(:,:) = zglam(:,:) - zlon 65 IF (jphgr_msh .NE. 2 .AND. jphgr_msh .NE. 3) THEN 66 zlon = MOD( plon + 720., 360. ) ! plon between 0 and 360 67 zglam(:,:) = MOD( zglam(:,:) + 720., 360. ) ! glam between 0 and 360 68 IF( zlon > 270. ) zlon = zlon - 360. ! zlon between -90 and 270 69 IF( zlon < 90. ) WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360. ! glam between -180 and 180 70 zglam(:,:) = zglam(:,:) - zlon 71 ELSE 72 zglam(:,:) = zglam(:,:) - plon 73 END IF 74 ! 67 75 zgphi(:,:) = zgphi(:,:) - plat 68 76 zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:) -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5506 r5619 28 28 USE in_out_manager ! I/O manager 29 29 USE iom ! I/O manager library 30 USE restart 30 USE restart , ONLY : rst_read_open ! ocean restart 31 31 USE lib_mpp ! distributed memory computing library 32 32 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 199 199 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 200 200 END DO 201 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )202 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )201 hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1. - ssumask(:,:) ) 202 hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1. - ssvmask(:,:) ) 203 203 204 204 ! Restoring frequencies for z_tilde coordinate … … 545 545 END DO 546 546 ! ! Inverse of the local depth 547 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)548 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)547 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 548 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 549 549 550 550 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5506 r5619 365 365 !! - bathy : meter bathymetry (in meters) 366 366 !!---------------------------------------------------------------------- 367 INTEGER :: ji, jj, j l, jk ! dummy loop indices367 INTEGER :: ji, jj, jk ! dummy loop indices 368 368 INTEGER :: inum ! temporary logical unit 369 369 INTEGER :: ierror ! error flag … … 472 472 risfdep(:,:)=0.e0 473 473 misfdep(:,:)=1 474 ! 475 ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 476 IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN 477 risfdep(:,:)=200.e0 478 misfdep(:,:)=1 479 ij0 = 1 ; ij1 = 40 480 DO jj = mj0(ij0), mj1(ij1) 481 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 482 END DO 483 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 484 ! 485 ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 486 ! 487 risfdep(:,:)=0.e0 488 misfdep(:,:)=1 489 ij0 = 1 ; ij1 = 40 490 DO jj = mj0(ij0), mj1(ij1) 491 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 492 END DO 493 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 494 END IF 474 495 ! 475 496 DEALLOCATE( idta, zdta ) … … 529 550 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 530 551 END IF 552 ! set grounded point to 0 553 WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 ) 554 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 555 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 556 END WHERE 531 557 ! 532 558 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration … … 952 978 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 953 979 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 954 REAL(wp) :: zmax ! Maximum depth955 980 REAL(wp) :: zdiff ! temporary scalar 956 REAL(wp) :: z refdep! temporary scalar981 REAL(wp) :: zmax ! temporary scalar 957 982 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 958 983 !!--------------------------------------------------------------------- … … 993 1018 END DO 994 1019 995 IF ( ln_isfcav ) CALL zgr_isf996 997 1020 ! Scale factors and depth at T- and W-points 998 1021 DO jk = 1, jpk ! intitialization to the reference z-coordinate … … 1002 1025 e3w_0 (:,:,jk) = e3w_1d (jk) 1003 1026 END DO 1004 ! 1005 DO jj = 1, jpj 1006 DO ji = 1, jpi 1007 ik = mbathy(ji,jj) 1008 IF( ik > 0 ) THEN ! ocean point only 1009 ! max ocean level case 1010 IF( ik == jpkm1 ) THEN 1011 zdepwp = bathy(ji,jj) 1012 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1013 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1014 e3t_0(ji,jj,ik ) = ze3tp 1015 e3t_0(ji,jj,ik+1) = ze3tp 1016 e3w_0(ji,jj,ik ) = ze3wp 1017 e3w_0(ji,jj,ik+1) = ze3tp 1018 gdepw_0(ji,jj,ik+1) = zdepwp 1019 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1020 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1021 ! 1022 ELSE ! standard case 1023 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1024 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1027 1028 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 1029 IF ( ln_isfcav ) CALL zgr_isf 1030 1031 ! Scale factors and depth at T- and W-points 1032 IF ( .NOT. ln_isfcav ) THEN 1033 DO jj = 1, jpj 1034 DO ji = 1, jpi 1035 ik = mbathy(ji,jj) 1036 IF( ik > 0 ) THEN ! ocean point only 1037 ! max ocean level case 1038 IF( ik == jpkm1 ) THEN 1039 zdepwp = bathy(ji,jj) 1040 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1041 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1042 e3t_0(ji,jj,ik ) = ze3tp 1043 e3t_0(ji,jj,ik+1) = ze3tp 1044 e3w_0(ji,jj,ik ) = ze3wp 1045 e3w_0(ji,jj,ik+1) = ze3tp 1046 gdepw_0(ji,jj,ik+1) = zdepwp 1047 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1048 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1049 ! 1050 ELSE ! standard case 1051 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1052 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1053 ENDIF 1054 !gm Bug? check the gdepw_1d 1055 ! ... on ik 1056 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1057 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1058 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1059 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1060 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1061 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1062 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1063 ! ... on ik+1 1064 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1065 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1066 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1025 1067 ENDIF 1026 !gm Bug? check the gdepw_1d1027 ! ... on ik1028 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &1029 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) &1030 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))1031 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &1032 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )1033 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &1034 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )1035 ! ... on ik+11036 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1037 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1038 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1039 1068 ENDIF 1040 END IF1041 END DO 1042 END DO1043 !1044 it = 01045 DO jj = 1, jpj1046 DO ji = 1, jpi1047 ik = mbathy(ji,jj)1048 IF( ik > 0 ) THEN ! ocean point only1049 e3tp (ji,jj) = e3t_0(ji,jj,ik)1050 e3wp (ji,jj) = e3w_0(ji,jj,ik)1051 ! test1052 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )1053 IF( zdiff <= 0._wp .AND. lwp ) THEN1054 it = it + 11055 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1056 WRITE(numout,*) ' bathy = ', bathy(ji,jj)1057 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1058 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )1069 END DO 1070 END DO 1071 ! 1072 it = 0 1073 DO jj = 1, jpj 1074 DO ji = 1, jpi 1075 ik = mbathy(ji,jj) 1076 IF( ik > 0 ) THEN ! ocean point only 1077 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1078 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1079 ! test 1080 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1081 IF( zdiff <= 0._wp .AND. lwp ) THEN 1082 it = it + 1 1083 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1084 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1085 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1086 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1087 ENDIF 1059 1088 ENDIF 1060 ENDIF 1061 END DO 1062 END DO 1063 ! 1064 IF ( ln_isfcav ) THEN 1065 ! (ISF) Definition of e3t, u, v, w for ISF case 1066 DO jj = 1, jpj 1067 DO ji = 1, jpi 1068 ik = misfdep(ji,jj) 1069 IF( ik > 1 ) THEN ! ice shelf point only 1070 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1071 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1072 !gm Bug? check the gdepw_0 1073 ! ... on ik 1074 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1075 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1076 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1077 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1078 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1079 1080 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1081 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1082 ENDIF 1083 ! ... on ik / ik-1 1084 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1085 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1086 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1087 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1088 ENDIF 1089 END DO 1090 END DO 1091 ! 1092 it = 0 1093 DO jj = 1, jpj 1094 DO ji = 1, jpi 1095 ik = misfdep(ji,jj) 1096 IF( ik > 1 ) THEN ! ice shelf point only 1097 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1098 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1099 ! test 1100 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1101 IF( zdiff <= 0. .AND. lwp ) THEN 1102 it = it + 1 1103 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1104 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1105 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1106 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1107 ENDIF 1108 ENDIF 1109 END DO 1110 END DO 1089 END DO 1090 END DO 1111 1091 END IF 1112 ! END (ISF) 1113 1092 ! 1114 1093 ! Scale factors and depth at U-, V-, UW and VW-points 1115 1094 DO jk = 1, jpk ! initialisation to z-scale factors … … 1119 1098 e3vw_0(:,:,jk) = e3w_1d(jk) 1120 1099 END DO 1100 1121 1101 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1122 1102 DO jj = 1, jpjm1 … … 1148 1128 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1149 1129 ! 1130 1150 1131 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1151 1132 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) … … 1255 1236 !!---------------------------------------------------------------------- 1256 1237 !! 1257 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1258 INTEGER :: ik, it ! temporary integers 1259 INTEGER :: id, jd, nprocd 1238 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1239 INTEGER :: ik, it ! temporary integers 1260 1240 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 1261 LOGICAL :: ll_print ! Allow control print for debugging 1241 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1242 REAL(wp) :: zmax ! Maximum and minimum depth 1243 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar 1262 1244 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1263 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 1264 REAL(wp) :: zmax, zmin ! Maximum and minimum depth 1245 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1265 1246 REAL(wp) :: zdiff ! temporary scalar 1266 REAL(wp) :: zrefdep ! temporary scalar1267 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar1268 1247 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1269 1248 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) … … 1280 1259 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1281 1260 END WHERE 1261 1262 ! set grounded point to 0 1263 WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 ) 1264 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1265 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1266 END WHERE 1282 1267 1283 1268 ! Compute misfdep for ocean points (i.e. first wet level) … … 1292 1277 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1293 1278 END WHERE 1279 1280 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1281 IF ( cp_cfg .NE. "isomip" ) THEN 1282 WHERE (risfdep(:,:) < 100 ) 1283 misfdep = 1; risfdep = 0.0_wp; 1284 END WHERE 1285 END IF 1294 1286 1295 1287 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation … … 1297 1289 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1298 1290 DO jl = 1, 10 1299 WHERE (bathy(:,:) . EQ. risfdep(:,:))1291 WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 ) 1300 1292 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1301 1293 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp … … 1326 1318 1327 1319 ! split last cell if possible (only where water column is 2 cell or less) 1328 DO jk = jpkm1, 1, -11329 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1330 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)1331 mbathy(:,:) = jk1332 bathy(:,:) = zmax1333 END WHERE1334 END DO1320 !DO jk = jpkm1, 1, -1 1321 ! zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1322 ! WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1323 ! mbathy(:,:) = jk 1324 ! bathy(:,:) = zmax 1325 ! END WHERE 1326 !END DO 1335 1327 1336 1328 ! split top cell if possible (only where water column is 2 cell or less) … … 1350 1342 ! test bathy 1351 1343 IF (risfdep(ji,jj) .GT. 1) THEN 1352 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1353 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1354 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1355 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1344 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1345 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1356 1346 1357 1347 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 1358 IF (zbathydiff .LE. zrisfdepdiff) THEN1359 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )1360 mbathy(ji,jj)= mbathy(ji,jj) + 11361 ELSE1348 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1349 ! bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1350 ! mbathy(ji,jj)= mbathy(ji,jj) + 1 1351 ! ELSE 1362 1352 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1363 1353 misfdep(ji,jj) = misfdep(ji,jj) - 1 1364 END IF1354 ! END IF 1365 1355 END IF 1366 1356 END IF … … 1374 1364 ! test bathy 1375 1365 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1376 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)& 1377 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1378 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1379 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1380 IF (zbathydiff .LE. zrisfdepdiff) THEN 1381 mbathy(ji,jj) = mbathy(ji,jj) + 1 1382 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1383 ELSE 1366 ! zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1367 ! zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1368 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1369 ! mbathy(ji,jj) = mbathy(ji,jj) + 1 1370 ! bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1371 ! ELSE 1384 1372 misfdep(ji,jj)= misfdep(ji,jj) - 1 1385 1373 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1386 END IF1374 ! END IF 1387 1375 ENDIF 1388 1376 END DO … … 1393 1381 DO ji = 1, jpim1 1394 1382 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1395 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) & 1396 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1397 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) & 1398 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1399 IF (zbathydiff .LE. zrisfdepdiff) THEN 1400 mbathy(ji,jj) = mbathy(ji,jj) + 1 1401 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) & 1402 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1403 ELSE 1383 ! zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1384 ! zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1385 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1386 ! mbathy(ji,jj) = mbathy(ji,jj) + 1 1387 ! bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1388 ! ELSE 1404 1389 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1405 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) & 1406 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1407 END IF 1390 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1391 ! END IF 1408 1392 ENDIF 1409 1393 END DO … … 1424 1408 DO ji = 1, jpim1 1425 1409 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 1426 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 1427 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1428 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) & 1429 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1430 IF (zbathydiff .LE. zrisfdepdiff) THEN 1431 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1432 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1433 ELSE 1410 ! zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1411 ! zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1412 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1413 ! mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1414 ! bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1415 ! ELSE 1434 1416 misfdep(ji,jj) = misfdep(ji,jj) - 1 1435 1417 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1436 END IF1418 ! END IF 1437 1419 ENDIF 1438 1420 END DO … … 1455 1437 DO ji = 1, jpim1 1456 1438 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1457 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1458 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1459 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 1460 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1461 IF (zbathydiff .LE. zrisfdepdiff) THEN 1462 mbathy(ji,jj) = mbathy(ji,jj) + 1 1463 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1464 ELSE 1439 ! zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1440 ! zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1441 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1442 ! mbathy(ji,jj) = mbathy(ji,jj) + 1 1443 ! bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1444 ! ELSE 1465 1445 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1466 1446 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1467 END IF1447 ! END IF 1468 1448 ENDIF 1469 1449 ENDDO … … 1485 1465 DO ji = 1, jpim1 1486 1466 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1487 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 1488 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1489 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) & 1490 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1491 IF (zbathydiff .LE. zrisfdepdiff) THEN 1492 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1493 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) & 1494 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1495 ELSE 1467 ! zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1468 ! zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1469 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1470 ! mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1471 ! bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1472 ! ELSE 1496 1473 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1497 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) & 1498 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1499 END IF 1474 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1475 ! END IF 1500 1476 ENDIF 1501 1477 ENDDO … … 1729 1705 ! end check compatibility ice shelf/bathy 1730 1706 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1731 WHERE (misfdep(:,:) <= 5)1732 misfdep = 1; risfdep = 0.0_wp;1733 END WHERE1734 1735 1707 IF( icompt == 0 ) THEN 1736 1708 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' … … 1738 1710 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1739 1711 ENDIF 1712 1713 ! compute scale factor and depth at T- and W- points 1714 DO jj = 1, jpj 1715 DO ji = 1, jpi 1716 ik = mbathy(ji,jj) 1717 IF( ik > 0 ) THEN ! ocean point only 1718 ! max ocean level case 1719 IF( ik == jpkm1 ) THEN 1720 zdepwp = bathy(ji,jj) 1721 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1722 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1723 e3t_0(ji,jj,ik ) = ze3tp 1724 e3t_0(ji,jj,ik+1) = ze3tp 1725 e3w_0(ji,jj,ik ) = ze3wp 1726 e3w_0(ji,jj,ik+1) = ze3tp 1727 gdepw_0(ji,jj,ik+1) = zdepwp 1728 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1729 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1730 ! 1731 ELSE ! standard case 1732 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1733 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1734 ENDIF 1735 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1736 !gm Bug? check the gdepw_1d 1737 ! ... on ik 1738 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1739 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1740 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1741 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1742 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1743 ! ... on ik+1 1744 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1745 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1746 ENDIF 1747 ENDIF 1748 END DO 1749 END DO 1750 ! 1751 it = 0 1752 DO jj = 1, jpj 1753 DO ji = 1, jpi 1754 ik = mbathy(ji,jj) 1755 IF( ik > 0 ) THEN ! ocean point only 1756 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1757 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1758 ! test 1759 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1760 IF( zdiff <= 0._wp .AND. lwp ) THEN 1761 it = it + 1 1762 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1763 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1764 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1765 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1766 ENDIF 1767 ENDIF 1768 END DO 1769 END DO 1770 ! 1771 ! (ISF) Definition of e3t, u, v, w for ISF case 1772 DO jj = 1, jpj 1773 DO ji = 1, jpi 1774 ik = misfdep(ji,jj) 1775 IF( ik > 1 ) THEN ! ice shelf point only 1776 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1777 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1778 !gm Bug? check the gdepw_0 1779 ! ... on ik 1780 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1781 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1782 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1783 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1784 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1785 1786 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1787 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1788 ENDIF 1789 ! ... on ik / ik-1 1790 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1791 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1792 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1793 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1794 ENDIF 1795 END DO 1796 END DO 1797 1798 it = 0 1799 DO jj = 1, jpj 1800 DO ji = 1, jpi 1801 ik = misfdep(ji,jj) 1802 IF( ik > 1 ) THEN ! ice shelf point only 1803 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1804 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1805 ! test 1806 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1807 IF( zdiff <= 0. .AND. lwp ) THEN 1808 it = it + 1 1809 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1810 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1811 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1812 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1813 ENDIF 1814 ENDIF 1815 END DO 1816 END DO 1740 1817 1741 1818 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r5215 r5619 34 34 35 35 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsd ! structure of input SST (file informations, fields read) 36 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsddmp ! structure of input SST (file informations, fields read) 36 37 37 38 !! * Substitutions … … 60 61 TYPE(FLD_N), DIMENSION( jpts) :: slf_i ! array of namelist informations on the fields to read 61 62 TYPE(FLD_N) :: sn_tem, sn_sal 63 TYPE(FLD_N) :: sn_dmpt, sn_dmps 62 64 !! 63 65 NAMELIST/namtsd/ ln_tsd_init, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal 66 NAMELIST/namtra_dmpfile/ sn_dmpt, sn_dmps 64 67 INTEGER :: ios 65 68 !!---------------------------------------------------------------------- … … 78 81 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in configuration namelist', lwp ) 79 82 IF(lwm) WRITE ( numond, namtsd ) 83 84 REWIND( numnam_ref ) ! Namelist namtra_dmp in reference namelist : Temperature and salinity damping term 85 READ ( numnam_ref, namtra_dmpfile, IOSTAT = ios, ERR = 903) 86 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp ) 87 88 REWIND( numnam_cfg ) ! Namelist namtra_dmp in configuration namelist : Temperature and salinity damping term 89 READ ( numnam_cfg, namtra_dmpfile, IOSTAT = ios, ERR = 904 ) 90 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp ) 80 91 81 92 IF( PRESENT( ld_tradmp ) ) ln_tsd_tradmp = .TRUE. ! forces the initialization when tradmp is used … … 105 116 ! 106 117 ALLOCATE( sf_tsd(jpts), STAT=ierr0 ) 118 ALLOCATE( sf_tsddmp(jpts), STAT=ierr0 ) 107 119 IF( ierr0 > 0 ) THEN 108 120 CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' ) ; RETURN … … 113 125 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk) , STAT=ierr2 ) 114 126 IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 127 ! dmp file 128 ALLOCATE( sf_tsddmp(jp_tem)%fnow(jpi,jpj,jpk) , STAT=ierr0 ) 129 IF( sn_dmpt%ln_tint ) ALLOCATE( sf_tsddmp(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 ) 130 ALLOCATE( sf_tsddmp(jp_sal)%fnow(jpi,jpj,jpk) , STAT=ierr2 ) 131 IF( sn_dmps%ln_tint ) ALLOCATE( sf_tsddmp(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 115 132 ! 116 133 IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN … … 120 137 slf_i(jp_tem) = sn_tem ; slf_i(jp_sal) = sn_sal 121 138 CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' ) 139 slf_i(jp_tem) = sn_dmpt ; slf_i(jp_sal) = sn_dmps 140 CALL fld_fill( sf_tsddmp, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' ) 122 141 ! 123 142 ENDIF … … 128 147 129 148 130 SUBROUTINE dta_tsd( kt, ptsd )149 SUBROUTINE dta_tsd( kt, ptsd, ptsddmp ) 131 150 !!---------------------------------------------------------------------- 132 151 !! *** ROUTINE dta_tsd *** … … 145 164 INTEGER , INTENT(in ) :: kt ! ocean time-step 146 165 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: ptsd ! T & S data 166 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), OPTIONAL, INTENT( out) :: ptsddmp ! T & S data 147 167 ! 148 168 INTEGER :: ji, jj, jk, jl, jkk ! dummy loop indicies … … 155 175 ! 156 176 CALL fld_read( kt, 1, sf_tsd ) !== read T & S data at kt time step ==! 177 IF ( PRESENT(ptsddmp) ) THEN 178 CALL fld_read( kt, 1, sf_tsddmp ) !== read T & S data at kt time step ==! 179 ptsddmp(:,:,:,jp_tem) = sf_tsddmp(jp_tem)%fnow(:,:,:) ! NO mask 180 ptsddmp(:,:,:,jp_sal) = sf_tsddmp(jp_sal)%fnow(:,:,:) 181 END IF 157 182 ! 158 183 ! … … 304 329 IF( sf_tsd(jp_sal)%ln_tint ) DEALLOCATE( sf_tsd(jp_sal)%fdta ) 305 330 DEALLOCATE( sf_tsd ) ! the structure itself 331 IF(lwp) WRITE(numout,*) 'dta_tsd: deallocte T & S arrays as they are only use to initialize the run' 332 DEALLOCATE( sf_tsddmp(jp_tem)%fnow ) ! T arrays in the structure 333 IF( sf_tsddmp(jp_tem)%ln_tint ) DEALLOCATE( sf_tsddmp(jp_tem)%fdta ) 334 DEALLOCATE( sf_tsddmp(jp_sal)%fnow ) ! S arrays in the structure 335 IF( sf_tsddmp(jp_sal)%ln_tint ) DEALLOCATE( sf_tsddmp(jp_sal)%fdta ) 336 DEALLOCATE( sf_tsddmp ) ! the structure itself 306 337 ENDIF 307 338 ! -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5332 r5619 47 47 USE wrk_nemo ! Memory allocation 48 48 USE timing ! Timing 49 USE sbc_iscpl 49 50 50 51 IMPLICIT NONE … … 91 92 ! ! ------------------- 92 93 CALL rst_read ! Read the restart file 94 IF (ln_iscpl) CALL rst_iscpl ! extraloate restart to wet and dry 93 95 CALL day_init ! model calendar (using both namelist and restart infos) 94 96 ELSE -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
r5147 r5619 73 73 REAL(wp), PUBLIC :: xlsn !: = lfus*rhosn (volumetric latent heat fusion of snow) [J/m3] 74 74 #else 75 REAL(wp), PUBLIC :: rhoic = 9 00._wp !: volumic mass of sea ice [kg/m3]75 REAL(wp), PUBLIC :: rhoic = 917._wp !: volumic mass of sea ice [kg/m3] 76 76 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: conductivity of the ice [W/m/K] 77 77 REAL(wp), PUBLIC :: rcpic = 1.8837e+6_wp !: volumetric specific heat for ice [J/m3/K]
Note: See TracChangeset
for help on using the changeset viewer.