Changeset 5619
- Timestamp:
- 2015-07-20T19:43:15+02:00 (9 years ago)
- Location:
- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM
- Files:
-
- 1 added
- 25 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/CONFIG/SHARED/field_def.xml
r5517 r5619 193 193 194 194 <!-- * variable related to ice shelf forcing * --> 195 <field id="fwfisf" long_name="Ice shelf melting" unit="kg/m2/s" /> 196 <field id="qisf" long_name="Ice Shelf Heat Flux" unit="W/m2" /> 197 <field id="isfgammat" long_name="transfert coefficient for isf (temperature)" unit="m/s" /> 198 <field id="isfgammas" long_name="transfert coefficient for isf (salinity)" unit="m/s" /> 199 <field id="stbl" long_name="salinity in the Losh tbl" unit="1e-3" /> 200 <field id="ttbl" long_name="temperature in the Losh tbl" unit="degC" /> 195 <field id="fwfisf" long_name="Ice shelf melting" unit="Kg/m2/s" /> 196 <field id="qisf" long_name="Ice Shelf Heat Flux" unit="W/m2" /> 197 <field id="isfgammat" long_name="transfert coefficient for isf (temperature) " unit="m/s" /> 198 <field id="isfgammas" long_name="transfert coefficient for isf (salinity) " unit="m/s" /> 199 <field id="stbl" long_name="salinity in the Losh tbl " unit="PSU" /> 200 <field id="ttbl" long_name="temperature in the Losh tbl " unit="C" /> 201 <field id="utbl" long_name="zonal current in the Losh tbl at T point " unit="m/s" /> 202 <field id="vtbl" long_name="merid current in the Losh tbl at T point " unit="m/s" /> 203 <field id="thermald" long_name="thermal driving of ice shelf melting " unit="C" /> 204 <field id="tfrz" long_name="top freezing point (used to compute melt) " unit="C" /> 205 <field id="tinsitu" long_name="top insitu temperature (used to cmpt melt) " unit="C" /> 206 <field id="ustar" long_name="ustar at T point used in ice shelf melting " unit="m/s" /> 201 207 202 208 <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/CONFIG/SHARED/namelist_ref
r5577 r5619 40 40 cn_ocerst_out = "restart" ! suffix of ocean restart name (output) 41 41 cn_ocerst_outdir = "." ! directory in which to write output ocean restarts 42 ln_iscpl = .false. ! cavity evolution forcing or coupling to ice sheet model 42 43 nn_istate = 0 ! output the initial state (1) or not (0) 43 44 ln_rst_list = .false. ! output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) … … 225 226 !! namsbc_rnf river runoffs 226 227 !! namsbc_isf ice shelf melting/freezing 228 !! namsbc_iscpl coupling option between land ice model and ocean 227 229 !! namsbc_apr Atmospheric Pressure 228 230 !! namsbc_ssr sea surface restoring term (for T and/or S) … … 469 471 / 470 472 !----------------------------------------------------------------------- 473 &namsbc_iscpl ! land ice / ocean coupling option 474 !----------------------------------------------------------------------- 475 rn_fiscpl = 43800 ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency) 476 ln_hfb = .false. ! activate conservation module (conservation exact after a time of rn_fiscpl) 477 / 478 !----------------------------------------------------------------------- 471 479 &namsbc_apr ! Atmospheric pressure used as ocean forcing or in bulk 472 480 !----------------------------------------------------------------------- … … 802 810 rn_smsh = 1. ! Smagorinsky diffusivity: = 0 - use only sheer 803 811 rn_aht_m = 2000. ! upper limit or stability criteria for lateral eddy diffusivity (m2/s) 812 / 813 !----------------------------------------------------------------------- 814 &namtra_dmpfile ! tracer: T & S newtonian damping 815 !----------------------------------------------------------------------- 816 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 817 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 818 sn_dmpt = 'resto', -1 ,'Tinit' , .true. , .true. , 'yearly' , '' , '' , '' 819 sn_dmps = 'resto', -1 ,'Sinit' , .true. , .true. , 'yearly' , '' , '' , '' 804 820 / 805 821 !----------------------------------------------------------------------- -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
r5586 r5619 197 197 ! Elevation 198 198 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*tmask_i(ji,jj) 199 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)* umask_i(ji,jj)200 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)* vmask_i(ji,jj)199 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 200 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 201 201 END DO 202 202 END DO … … 326 326 X1= ana_amp(ji,jj,jh,1) 327 327 X2=-ana_amp(ji,jj,jh,2) 328 out_u(ji,jj, jh) = X1 * umask_i(ji,jj)329 out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)328 out_u(ji,jj, jh) = X1 * ssumask(ji,jj) 329 out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj) 330 330 ENDDO 331 331 ENDDO -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
r5120 r5619 165 165 DO jk = 1, jpkm1 166 166 ! volume variation (calculated with scale factors) 167 zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk) &168 & * (fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )167 zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * ( tmask(:,:,jk) & 168 & * fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) ) 169 169 ! heat content variation 170 zdiff_hc = zdiff_hc + glob_sum( surf(:,:) * tmask(:,:,jk) &171 & * (fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) )170 zdiff_hc = zdiff_hc + glob_sum( surf(:,:) * ( tmask(:,:,jk) & 171 & * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) ) 172 172 ! salt content variation 173 zdiff_sc = zdiff_sc + glob_sum( surf(:,:) * tmask(:,:,jk) &174 & * (fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) )173 zdiff_sc = zdiff_sc + glob_sum( surf(:,:) * ( tmask(:,:,jk) & 174 & * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) ) 175 175 ENDDO 176 176 … … 200 200 ! ENDIF 201 201 !!gm end 202 202 203 IF (lwp) PRINT *, 'ISCPL CONS HEAT ', kt, zdiff_hc / zvol_tot, zdiff_sc / zvol_tot 204 IF (lwp) PRINT *, 'ISCPL CONS VOL ', kt, zdiff_v1 * 1.e-9, zdiff_v2 * 1.e-9 203 205 204 206 IF( lk_vvl ) THEN … … 217 219 CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp ) ! Heat content variation (1.e20 J) 218 220 CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9 ) ! Salt content variation (psu*km3) 219 CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 )! volume ssh variation (km3)221 CALL iom_put( 'bgvolssh' , (zdiff_v1+zdiff_v2) * 1.e-9 ) ! volume ssh variation (km3) 220 222 CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 ) ! vol - surface forcing (km3) 221 223 CALL iom_put( 'bgfrctem' , frc_t / zvol_tot ) ! hc - surface forcing (C) … … 277 279 ssh_ini(:,:) = sshn(:,:) ! initial ssh 278 280 DO jk = 1, jpk 279 e3t_ini (:,:,jk) = fse3t_n(:,:,jk) ! initial vertical scale factors280 hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) ! initial heat content281 sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) ! initial salt content281 e3t_ini (:,:,jk) = fse3t_n(:,:,jk) * tmask(:,:,jk) ! initial vertical scale factors 282 hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) * tmask(:,:,jk) ! initial heat content 283 sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) * tmask(:,:,jk) ! initial salt content 282 284 END DO 283 285 frc_v = 0._wp ! volume trend due to forcing -
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] -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
r5516 r5619 29 29 USE sbcrnf ! river runoff 30 30 USE sbcisf ! ice shelf 31 USE sbc_iscpl ! ice shelf 31 32 USE cla ! cross land advection (cla_div routine) 32 33 USE in_out_manager ! I/O manager … … 328 329 329 330 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 330 IF( ln_divisf .AND. (nn_isf .GT. 0) ) CALL sbc_isf_div( hdivn ) ! ice shelf (update hdivn field) 331 IF( ln_divisf .AND. (nn_isf .GT. 0) ) CALL sbc_isf_div ( hdivn ) ! ice shelf (update hdivn field) 332 IF( ln_iscpl .AND. ln_hfb ) CALL sbc_iscpl_div( hdivn ) ! ice shelf (update hdivn field) 331 333 IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 332 334 ! -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5467 r5619 353 353 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 354 354 END DO 355 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )356 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )355 hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 356 hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 357 357 ENDIF 358 358 ! -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
r5407 r5619 24 24 USE trdmxl_oce ! ocean active mixed layer tracers trends variables 25 25 USE divcur ! hor. divergence and curl (div & cur routines) 26 USE wrk_nemo 26 27 27 28 IMPLICIT NONE … … 145 146 CALL iom_rstput( kt, nitrst, numrow, 'rhd' , rhd ) 146 147 #endif 148 149 150 151 IF ( ln_iscpl ) THEN 152 CALL iom_rstput( kt, nitrst, numrow, 'tmask' , tmask ) ! need to extrapolate T/S 153 CALL iom_rstput( kt, nitrst, numrow, 'umask' , umask ) ! need to correct barotropic velocity 154 CALL iom_rstput( kt, nitrst, numrow, 'vmask' , vmask ) ! need to correct barotropic velocity 155 CALL iom_rstput( kt, nitrst, numrow, 'smask' , ssmask ) ! need to correct barotropic velocity 156 CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) ) ! need to compute temperature correction 157 CALL iom_rstput( kt, nitrst, numrow, 'fse3u_n', fse3u_n(:,:,:) ) ! need to compute volume correction ???? 158 CALL iom_rstput( kt, nitrst, numrow, 'fse3v_n', fse3v_n(:,:,:) ) ! need to compute volume correction ???? 159 CALL iom_rstput( kt, nitrst, numrow, 'fsdepw_n', fsdepw_n(:,:,:) ) ! need to compute volume correction ???? 160 END IF 147 161 IF( kt == nitrst ) THEN 148 162 CALL iom_close( numrow ) ! close the restart file (only at last time step) -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90
r5429 r5619 31 31 END INTERFACE 32 32 33 INTERFACE lbc_sum 34 MODULE PROCEDURE mpp_sum_3d, mpp_sum_2d 35 END INTERFACE 36 33 37 INTERFACE lbc_bdy_lnk 34 38 MODULE PROCEDURE mpp_lnk_bdy_2d, mpp_lnk_bdy_3d … … 45 49 PUBLIC lbc_lnk ! ocean lateral boundary conditions 46 50 PUBLIC lbc_lnk_multi ! modified ocean lateral boundary conditions 51 PUBLIC lbc_sum 47 52 PUBLIC lbc_lnk_e 48 53 PUBLIC lbc_bdy_lnk ! ocean lateral BDY boundary conditions -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
r5429 r5619 72 72 PUBLIC mpp_lnk_3d, mpp_lnk_3d_gather, mpp_lnk_2d, mpp_lnk_2d_e 73 73 PUBLIC mpp_lnk_2d_9 74 PUBLIC mpp_sum_3d, mpp_sum_2d 74 75 PUBLIC mppscatter, mppgather 75 76 PUBLIC mpp_ini_ice, mpp_ini_znl … … 1394 1395 END SUBROUTINE mpp_lnk_2d_e 1395 1396 1397 SUBROUTINE mpp_sum_3d( ptab, cd_type, psgn, cd_mpp, pval ) 1398 !!---------------------------------------------------------------------- 1399 !! *** routine mpp_sum_3d *** 1400 !! 1401 !! ** Purpose : Message passing manadgement 1402 !! 1403 !! ** Method : Use mppsend and mpprecv function for passing mask 1404 !! between processors following neighboring subdomains. 1405 !! domain parameters 1406 !! nlci : first dimension of the local subdomain 1407 !! nlcj : second dimension of the local subdomain 1408 !! nbondi : mark for "east-west local boundary" 1409 !! nbondj : mark for "north-south local boundary" 1410 !! noea : number for local neighboring processors 1411 !! nowe : number for local neighboring processors 1412 !! noso : number for local neighboring processors 1413 !! nono : number for local neighboring processors 1414 !! 1415 !! ** Action : ptab with update value at its periphery 1416 !! 1417 !!---------------------------------------------------------------------- 1418 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: ptab ! 3D array on which the boundary condition is applied 1419 CHARACTER(len=1) , INTENT(in ) :: cd_type ! define the nature of ptab array grid-points 1420 ! ! = T , U , V , F , W points 1421 REAL(wp) , INTENT(in ) :: psgn ! =-1 the sign change across the north fold boundary 1422 ! ! = 1. , the sign is kept 1423 CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only 1424 REAL(wp) , OPTIONAL , INTENT(in ) :: pval ! background value (used at closed boundaries) 1425 !! 1426 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1427 INTEGER :: imigr, iihom, ijhom ! temporary integers 1428 INTEGER :: ml_req1, ml_req2, ml_err ! for key_mpi_isend 1429 REAL(wp) :: zland 1430 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for key_mpi_isend 1431 ! 1432 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zt3ns, zt3sn ! 3d for north-south & south-north 1433 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zt3ew, zt3we ! 3d for east-west & west-east 1434 1435 !!---------------------------------------------------------------------- 1436 1437 ALLOCATE( zt3ns(jpi,jprecj,jpk,2), zt3sn(jpi,jprecj,jpk,2), & 1438 & zt3ew(jpj,jpreci,jpk,2), zt3we(jpj,jpreci,jpk,2) ) 1439 1440 ! 1441 IF( PRESENT( pval ) ) THEN ; zland = pval ! set land value 1442 ELSE ; zland = 0.e0 ! zero by default 1443 ENDIF 1444 1445 ! 1. standard boundary treatment 1446 ! ------------------------------ 1447 IF( PRESENT( cd_mpp ) ) THEN ! only fill added line/raw with existing values 1448 ! 1449 ! WARNING ptab is defined only between nld and nle 1450 ! DO jk = 1, jpk 1451 ! DO jj = nlcj+1, jpj ! added line(s) (inner only) 1452 ! ptab(nldi :nlei , jj ,jk) = ptab(nldi:nlei, nlej,jk) 1453 ! ptab(1 :nldi-1, jj ,jk) = ptab(nldi , nlej,jk) 1454 ! ptab(nlei+1:nlci , jj ,jk) = ptab( nlei, nlej,jk) 1455 ! END DO 1456 ! DO ji = nlci+1, jpi ! added column(s) (full) 1457 ! ptab(ji ,nldj :nlej ,jk) = ptab( nlei,nldj:nlej,jk) 1458 ! ptab(ji ,1 :nldj-1,jk) = ptab( nlei,nldj ,jk) 1459 ! ptab(ji ,nlej+1:jpj ,jk) = ptab( nlei, nlej,jk) 1460 ! END DO 1461 ! END DO 1462 ! 1463 ELSE ! standard close or cyclic treatment 1464 ! 1465 ! ! East-West boundaries 1466 ! !* Cyclic east-west 1467 IF( nbondi == 2 .AND. (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN 1468 ! ptab( 1 ,:,:) = ptab(jpim1,:,:) 1469 ! ptab(jpi,:,:) = ptab( 2 ,:,:) 1470 ELSE !* closed 1471 ! IF( .NOT. cd_type == 'F' ) ptab( 1 :jpreci,:,:) = zland ! south except F-point 1472 ! ptab(nlci-jpreci+1:jpi ,:,:) = zland ! north 1473 ENDIF 1474 ! ! North-South boundaries (always closed) 1475 ! IF( .NOT. cd_type == 'F' ) ptab(:, 1 :jprecj,:) = zland ! south except F-point 1476 ! ptab(:,nlcj-jprecj+1:jpj ,:) = zland ! north 1477 ! 1478 ENDIF 1479 1480 ! 2. East and west directions exchange 1481 ! ------------------------------------ 1482 ! we play with the neigbours AND the row number because of the periodicity 1483 ! 1484 SELECT CASE ( nbondi ) ! Read Dirichlet lateral conditions 1485 CASE ( -1, 0, 1 ) ! all exept 2 (i.e. close case) 1486 iihom = nlci-jpreci 1487 DO jl = 1, jpreci 1488 zt3ew(:,jl,:,1) = ptab(jl ,:,:) ; ptab(jl ,:,:) = 0.0_wp 1489 zt3we(:,jl,:,1) = ptab(iihom+jl,:,:) ; ptab(iihom+jl,:,:) = 0.0_wp 1490 END DO 1491 END SELECT 1492 ! 1493 ! ! Migrations 1494 imigr = jpreci * jpj * jpk 1495 ! 1496 SELECT CASE ( nbondi ) 1497 CASE ( -1 ) 1498 CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req1 ) 1499 CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea ) 1500 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1501 CASE ( 0 ) 1502 CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 ) 1503 CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req2 ) 1504 CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea ) 1505 CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe ) 1506 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1507 IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 1508 CASE ( 1 ) 1509 CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 ) 1510 CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe ) 1511 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1512 END SELECT 1513 ! 1514 ! ! Write Dirichlet lateral conditions 1515 iihom = nlci-nreci 1516 ! 1517 SELECT CASE ( nbondi ) 1518 CASE ( -1 ) 1519 DO jl = 1, jpreci 1520 ptab(iihom+jl,:,:) = ptab(iihom+jl,:,:) + zt3ew(:,jl,:,2) 1521 END DO 1522 CASE ( 0 ) 1523 DO jl = 1, jpreci 1524 ptab(jpreci+jl,:,:) = ptab(jpreci+jl,:,:) + zt3we(:,jl,:,2) 1525 ptab(iihom +jl,:,:) = ptab(iihom +jl,:,:) + zt3ew(:,jl,:,2) 1526 END DO 1527 CASE ( 1 ) 1528 DO jl = 1, jpreci 1529 ptab(jpreci+jl,:,:) = ptab(jpreci+jl,:,:) + zt3we(:,jl,:,2) 1530 END DO 1531 END SELECT 1532 1533 1534 ! 3. North and south directions 1535 ! ----------------------------- 1536 ! always closed : we play only with the neigbours 1537 ! 1538 IF( nbondj /= 2 ) THEN ! Read Dirichlet lateral conditions 1539 ijhom = nlcj-jprecj 1540 DO jl = 1, jprecj 1541 zt3sn(:,jl,:,1) = ptab(:,ijhom+jl,:) ; ptab(:,ijhom+jl,:) = 0.0_wp 1542 zt3ns(:,jl,:,1) = ptab(:,jl ,:) ; ptab(:,jl ,:) = 0.0_wp 1543 END DO 1544 ENDIF 1545 ! 1546 ! ! Migrations 1547 imigr = jprecj * jpi * jpk 1548 ! 1549 SELECT CASE ( nbondj ) 1550 CASE ( -1 ) 1551 CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req1 ) 1552 CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono ) 1553 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1554 CASE ( 0 ) 1555 CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 ) 1556 CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req2 ) 1557 CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono ) 1558 CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso ) 1559 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1560 IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 1561 CASE ( 1 ) 1562 CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 ) 1563 CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso ) 1564 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 1565 END SELECT 1566 ! 1567 ! ! Write Dirichlet lateral conditions 1568 ijhom = nlcj-nrecj 1569 ! 1570 SELECT CASE ( nbondj ) 1571 CASE ( -1 ) 1572 DO jl = 1, jprecj 1573 ptab(:,ijhom+jl,:) = ptab(:,ijhom+jl,:) + zt3ns(:,jl,:,2) 1574 END DO 1575 CASE ( 0 ) 1576 DO jl = 1, jprecj 1577 ptab(:,jprecj+jl,:) = ptab(:,jprecj+jl,:) + zt3sn(:,jl,:,2) 1578 ptab(:,ijhom +jl,:) = ptab(:,ijhom +jl,:) + zt3ns(:,jl,:,2) 1579 END DO 1580 CASE ( 1 ) 1581 DO jl = 1, jprecj 1582 ptab(:,jprecj+jl,:) = ptab(:,jprecj+jl,:) + zt3sn(:,jl ,:,2) 1583 END DO 1584 END SELECT 1585 1586 1587 ! 4. north fold treatment 1588 ! ----------------------- 1589 ! 1590 IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN 1591 ! 1592 SELECT CASE ( jpni ) 1593 CASE ( 1 ) ; CALL lbc_nfd ( ptab, cd_type, psgn ) ! only 1 northern proc, no mpp 1594 CASE DEFAULT ; CALL mpp_lbc_north( ptab, cd_type, psgn ) ! for all northern procs. 1595 END SELECT 1596 ! 1597 ENDIF 1598 ! 1599 DEALLOCATE( zt3ns, zt3sn, zt3ew, zt3we ) 1600 ! 1601 END SUBROUTINE mpp_sum_3d 1602 1603 SUBROUTINE mpp_sum_2d( pt2d, cd_type, psgn, cd_mpp, pval ) 1604 !!---------------------------------------------------------------------- 1605 !! *** routine mpp_sum_2d *** 1606 !! 1607 !! ** Purpose : Message passing manadgement for 2d array 1608 !! 1609 !! ** Method : Use mppsend and mpprecv function for passing mask 1610 !! between processors following neighboring subdomains. 1611 !! domain parameters 1612 !! nlci : first dimension of the local subdomain 1613 !! nlcj : second dimension of the local subdomain 1614 !! nbondi : mark for "east-west local boundary" 1615 !! nbondj : mark for "north-south local boundary" 1616 !! noea : number for local neighboring processors 1617 !! nowe : number for local neighboring processors 1618 !! noso : number for local neighboring processors 1619 !! nono : number for local neighboring processors 1620 !! 1621 !!---------------------------------------------------------------------- 1622 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pt2d ! 2D array on which the boundary condition is applied 1623 CHARACTER(len=1) , INTENT(in ) :: cd_type ! define the nature of ptab array grid-points 1624 ! ! = T , U , V , F , W and I points 1625 REAL(wp) , INTENT(in ) :: psgn ! =-1 the sign change across the north fold boundary 1626 ! ! = 1. , the sign is kept 1627 CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only 1628 REAL(wp) , OPTIONAL , INTENT(in ) :: pval ! background value (used at closed boundaries) 1629 !! 1630 INTEGER :: ji, jj, jl ! dummy loop indices 1631 INTEGER :: imigr, iihom, ijhom ! temporary integers 1632 INTEGER :: ml_req1, ml_req2, ml_err ! for key_mpi_isend 1633 REAL(wp) :: zland 1634 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for key_mpi_isend 1635 ! 1636 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zt2ns, zt2sn ! 2d for north-south & south-north 1637 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zt2ew, zt2we ! 2d for east-west & west-east 1638 1639 !!---------------------------------------------------------------------- 1640 1641 ALLOCATE( zt2ns(jpi,jprecj,2), zt2sn(jpi,jprecj,2), & 1642 & zt2ew(jpj,jpreci,2), zt2we(jpj,jpreci,2) ) 1643 1644 ! 1645 IF( PRESENT( pval ) ) THEN ; zland = pval ! set land value 1646 ELSE ; zland = 0.e0 ! zero by default 1647 ENDIF 1648 1649 ! 1. standard boundary treatment 1650 ! ------------------------------ 1651 ! 1652 ! IF( PRESENT( cd_mpp ) ) THEN ! only fill added line/raw with existing values 1653 ! ! 1654 ! ! WARNING pt2d is defined only between nld and nle 1655 ! DO jj = nlcj+1, jpj ! added line(s) (inner only) 1656 ! pt2d(nldi :nlei , jj ) = pt2d(nldi:nlei, nlej) 1657 ! pt2d(1 :nldi-1, jj ) = pt2d(nldi , nlej) 1658 ! pt2d(nlei+1:nlci , jj ) = pt2d( nlei, nlej) 1659 ! END DO 1660 ! DO ji = nlci+1, jpi ! added column(s) (full) 1661 ! pt2d(ji ,nldj :nlej ) = pt2d( nlei,nldj:nlej) 1662 ! pt2d(ji ,1 :nldj-1) = pt2d( nlei,nldj ) 1663 ! pt2d(ji ,nlej+1:jpj ) = pt2d( nlei, nlej) 1664 ! END DO 1665 ! ! 1666 ! ELSE ! standard close or cyclic treatment 1667 ! ! 1668 ! ! ! East-West boundaries 1669 ! IF( nbondi == 2 .AND. & ! Cyclic east-west 1670 ! & (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN 1671 ! pt2d( 1 ,:) = pt2d(jpim1,:) ! west 1672 ! pt2d(jpi,:) = pt2d( 2 ,:) ! east 1673 ! ELSE ! closed 1674 ! IF( .NOT. cd_type == 'F' ) pt2d( 1 :jpreci,:) = zland ! south except F-point 1675 ! pt2d(nlci-jpreci+1:jpi ,:) = zland ! north 1676 ! ENDIF 1677 ! ! ! North-South boundaries (always closed) 1678 ! IF( .NOT. cd_type == 'F' ) pt2d(:, 1 :jprecj) = zland !south except F-point 1679 ! pt2d(:,nlcj-jprecj+1:jpj ) = zland ! north 1680 ! ! 1681 ! ENDIF 1682 1683 ! 2. East and west directions exchange 1684 ! ------------------------------------ 1685 ! we play with the neigbours AND the row number because of the periodicity 1686 ! 1687 SELECT CASE ( nbondi ) ! Read Dirichlet lateral conditions 1688 CASE ( -1, 0, 1 ) ! all exept 2 (i.e. close case) 1689 iihom = nlci - jpreci 1690 DO jl = 1, jpreci 1691 zt2ew(:,jl,1) = pt2d(jl ,:) ; pt2d(jl ,:) = 0.0_wp 1692 zt2we(:,jl,1) = pt2d(iihom +jl,:) ; pt2d(iihom +jl,:) = 0.0_wp 1693 END DO 1694 END SELECT 1695 ! 1696 ! ! Migrations 1697 imigr = jpreci * jpj 1698 ! 1699 SELECT CASE ( nbondi ) 1700 CASE ( -1 ) 1701 CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req1 ) 1702 CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea ) 1703 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1704 CASE ( 0 ) 1705 CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 ) 1706 CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req2 ) 1707 CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea ) 1708 CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe ) 1709 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1710 IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 1711 CASE ( 1 ) 1712 CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 ) 1713 CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe ) 1714 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1715 END SELECT 1716 ! 1717 ! ! Write Dirichlet lateral conditions 1718 iihom = nlci-nreci 1719 ! 1720 SELECT CASE ( nbondi ) 1721 CASE ( -1 ) 1722 DO jl = 1, jpreci 1723 pt2d(iihom+jl,:) = pt2d(iihom+jl,:) + zt2ew(:,jl,2) 1724 END DO 1725 CASE ( 0 ) 1726 DO jl = 1, jpreci 1727 pt2d(jpreci+jl,:) = pt2d(jpreci+jl,:) + zt2we(:,jl,2) 1728 pt2d(iihom +jl,:) = pt2d(iihom +jl,:) + zt2ew(:,jl,2) 1729 END DO 1730 CASE ( 1 ) 1731 DO jl = 1, jpreci 1732 pt2d(jpreci+jl,:) = pt2d(jpreci+jl,:) + zt2we(:,jl,2) 1733 END DO 1734 END SELECT 1735 1736 1737 ! 3. North and south directions 1738 ! ----------------------------- 1739 ! always closed : we play only with the neigbours 1740 ! 1741 IF( nbondj /= 2 ) THEN ! Read Dirichlet lateral conditions 1742 ijhom = nlcj - jprecj 1743 DO jl = 1, jprecj 1744 zt2sn(:,jl,1) = pt2d(:,ijhom +jl) ; pt2d(:,ijhom +jl) = 0.0_wp 1745 zt2ns(:,jl,1) = pt2d(:,jl ) ; pt2d(:,jl ) = 0.0_wp 1746 END DO 1747 ENDIF 1748 ! 1749 ! ! Migrations 1750 imigr = jprecj * jpi 1751 ! 1752 SELECT CASE ( nbondj ) 1753 CASE ( -1 ) 1754 CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req1 ) 1755 CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono ) 1756 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1757 CASE ( 0 ) 1758 CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 ) 1759 CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req2 ) 1760 CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono ) 1761 CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso ) 1762 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1763 IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 1764 CASE ( 1 ) 1765 CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 ) 1766 CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso ) 1767 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1768 END SELECT 1769 ! 1770 ! ! Write Dirichlet lateral conditions 1771 ijhom = nlcj-nrecj 1772 ! 1773 SELECT CASE ( nbondj ) 1774 CASE ( -1 ) 1775 DO jl = 1, jprecj 1776 pt2d(:,ijhom+jl) = pt2d(:,ijhom+jl) + zt2ns(:,jl,2) 1777 END DO 1778 CASE ( 0 ) 1779 DO jl = 1, jprecj 1780 pt2d(:,jprecj+jl) = pt2d(:,jprecj+jl) + zt2sn(:,jl,2) 1781 pt2d(:,ijhom +jl) = pt2d(:,ijhom +jl) + zt2ns(:,jl,2) 1782 END DO 1783 CASE ( 1 ) 1784 DO jl = 1, jprecj 1785 pt2d(:,jprecj+jl) = pt2d(:,jprecj+jl) + zt2sn(:,jl,2) 1786 END DO 1787 END SELECT 1788 1789 1790 ! 4. north fold treatment 1791 ! ----------------------- 1792 ! 1793 IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN 1794 ! 1795 SELECT CASE ( jpni ) 1796 CASE ( 1 ) ; CALL lbc_nfd ( pt2d, cd_type, psgn ) ! only 1 northern proc, no mpp 1797 CASE DEFAULT ; CALL mpp_lbc_north( pt2d, cd_type, psgn ) ! for all northern procs. 1798 END SELECT 1799 ! 1800 ENDIF 1801 ! 1802 DEALLOCATE( zt2ns, zt2sn, zt2ew, zt2we ) 1803 ! 1804 END SUBROUTINE mpp_sum_2d 1396 1805 1397 1806 SUBROUTINE mppsend( ktyp, pmess, kbytes, kdest, md_req ) -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90
r5130 r5619 136 136 137 137 imask(:,:)=1 138 WHERE ( zdta(:,:) - zdtaisf(:,:) <= 0.) imask = 0138 WHERE ( zdta(:,:) - zdtaisf(:,:) <= 1e-2 ) imask = 0 139 139 140 140 ! 1. Dimension arrays for subdomains -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
r5120 r5619 111 111 zcoef = z_fwf * rcp 112 112 emp(:,:) = emp(:,:) - z_fwf * tmask(:,:,1) 113 sfx(:,:) = sfx(:,:) + z_fwf * sss_m * tmask(:,:,1) 113 114 qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction 114 115 ENDIF -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r5541 r5619 63 63 64 64 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! phycst ? 65 REAL(wp), PUBLIC, SAVE :: kappa = 1.54e-6_wp ! phycst ?65 REAL(wp), PUBLIC, SAVE :: kappa = 0.0_wp ! phycst ? 66 66 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! phycst ? 67 67 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! phycst ? … … 125 125 IF ( lwp ) WRITE(numout,*) ' ln_divisf = ', ln_divisf 126 126 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 127 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 128 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 127 129 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 128 130 IF (ln_divisf) THEN ! keep it in the namelist ??? used true anyway as for runoff ? (PM) … … 150 152 !: read effective lenght (BG03) 151 153 IF (nn_isf == 2) THEN 152 ! Read Data and save some integral values154 cvarLeff = 'soLeff' 153 155 CALL iom_open( sn_Leff_isf%clname, inum ) 154 cvarLeff = 'soLeff' !: variable name for Efficient Length scale155 156 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 156 157 CALL iom_close(inum) … … 237 238 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 238 239 CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 239 CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U')240 CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V')240 CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U') 241 CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V') 241 242 ! iom print 242 243 CALL iom_put('ttbl',ttbl(:,:)) 243 244 CALL iom_put('stbl',stbl(:,:)) 244 CALL iom_put('utbl',utbl(:,:) )245 CALL iom_put('vtbl',vtbl(:,:) )245 CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 246 CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 246 247 ! compute fwf and heat flux 247 248 CALL sbc_isf_cav (kt) … … 296 297 ! 297 298 ! output 298 CALL iom_put('qisf' , qisf)299 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce)299 IF( iom_use('qisf' ) ) CALL iom_put('qisf' , qisf) 300 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf) 300 301 END IF 301 302 … … 416 417 REAL(wp) :: zfwflx, zhtflx, zhtflx_b 417 418 REAL(wp) :: zgammat, zgammas 418 REAL(wp) :: zeps = -1.e-20_wp !== Local constant initialization ==!419 REAL(wp) :: zeps = 1.e-20_wp !== Local constant initialization ==! 419 420 INTEGER :: ji, jj ! dummy loop indices 420 421 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers … … 506 507 zeps1=rcp*rau0*zgammat 507 508 zeps2=lfusisf*rau0*zgammas 508 zeps3=rhoisf*rcpi*kappa/ risfdep(ji,jj)509 zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps) 509 510 zeps4=zlamb2+zlamb3*risfdep(ji,jj) 510 511 zeps6=zeps4-zti(ji,jj) 511 512 zeps7=zeps4-tsurf 512 513 zaqe=zlamb1 * (zeps1 + zeps3) 513 zaqer=0.5/ zaqe514 zaqer=0.5/MIN(zaqe,-zeps) 514 515 zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 515 516 zcqe=zeps2*stbl(ji,jj) … … 520 521 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 521 522 522 ! zfwflx is upward water flux 523 zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz)523 ! zfwflx is upward water flux (MAX usefull if kappa = 0 524 zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) ) 524 525 ! zhtflx is upward heat flux (out of ocean) 525 526 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value … … 527 528 ! zwflx is upward water flux 528 529 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 529 zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)530 !!!!!!!! zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj) 530 531 ! test convergence and compute gammat 531 532 IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r5541 r5619 921 921 pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 922 922 & - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) & 923 & / fse3w(ji,jj,jk) * tmask(ji,jj,jk)923 & / fse3w(ji,jj,jk) * wmask(ji,jj,jk) 924 924 END DO 925 925 END DO … … 1242 1242 ! 1243 1243 rau0 = 1026._wp !: volumic mass of reference [kg/m3] 1244 rcp = 39 91.86795711963_wp !: heat capacity [J/K]1244 rcp = 3974._wp !: heat capacity [J/K] 1245 1245 ! 1246 1246 IF(lwp) THEN ! Control print -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
r5102 r5619 102 102 REAL(wp) :: zta, zsa ! local scalars 103 103 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zts_dta 104 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zts_dtadmp 104 105 !!---------------------------------------------------------------------- 105 106 ! 106 107 IF( nn_timing == 1 ) CALL timing_start( 'tra_dmp') 107 108 ! 108 CALL wrk_alloc( jpi, jpj, jpk, jpts, zts_dta ) 109 ! 109 CALL wrk_alloc( jpi, jpj, jpk, jpts, zts_dta, zts_dtadmp ) 110 110 ! !== input T-S data at kt ==! 111 CALL dta_tsd( kt, zts_dta ) ! read and interpolates T-S data at kt 111 CALL dta_tsd( kt, zts_dta, zts_dtadmp ) ! read and interpolates T-S data at kt 112 zts_dta=zts_dtadmp 112 113 ! 113 114 SELECT CASE ( nn_zdmp ) !== type of damping ==! … … 175 176 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 176 177 ! 177 CALL wrk_dealloc( jpi, jpj, jpk, jpts, zts_dta )178 CALL wrk_dealloc( jpi, jpj, jpk, jpts, zts_dta, zts_dtadmp ) 178 179 ! 179 180 IF( nn_timing == 1 ) CALL timing_stop( 'tra_dmp') -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
r5431 r5619 22 22 USE sbcrnf ! River runoff 23 23 USE sbcisf ! Ice shelf 24 USE sbc_iscpl ! Ice sheet coupling 24 25 USE traqsr ! solar radiation penetration 25 26 USE trd_oce ! trends: ocean variables … … 117 118 INTEGER :: ji, jj, jk, jn ! dummy loop indices 118 119 INTEGER :: ikt, ikb 119 INTEGER :: nk_isf120 120 REAL(wp) :: zfact, z1_e3t, zdep 121 REAL(wp) :: zalpha, zhk 122 REAL(wp) :: zt_frz, zpress 121 REAL(wp) :: zt_frz, zpress 123 122 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds 124 123 !!---------------------------------------------------------------------- … … 220 219 ! 221 220 IF( nn_isf > 0 ) THEN 222 zfact = 0.5 e0221 zfact = 0.5_wp 223 222 DO jj = 2, jpj 224 223 DO ji = fs_2, fs_jpim1 … … 233 232 ! compute tfreez for the temperature correction (we add water at freezing temperature) 234 233 ! zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 235 zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )234 zt_frz = -1.9_wp !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 236 235 ! compute trend 237 236 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & … … 246 245 ! compute tfreez for the temperature correction (we add water at freezing temperature) 247 246 ! zpress = grav*rau0*fsdept(ji,jj,ikb)*1.e-04 248 zt_frz = -1.9 !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress )247 zt_frz = -1.9_wp !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress ) 249 248 ! compute trend 250 249 tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem) & … … 288 287 ENDIF 289 288 290 IF( l_trdtra ) THEN ! send trends for further diagnostics 289 !---------------------------------------- 290 ! Ice Sheet coupling imbalance correction to have conservation 291 !---------------------------------------- 292 ! 293 IF( ln_iscpl .AND. ln_hfb) THEN ! input of heat and salt due to river runoff 294 DO jk = 1,jpk 295 DO jj = 2, jpj 296 DO ji = fs_2, fs_jpim1 297 zdep = 1._wp / fse3t_n(ji,jj,jk) 298 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) & 299 & * zdep 300 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) & 301 & * zdep 302 END DO 303 END DO 304 END DO 305 ENDIF 306 307 IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics 291 308 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 292 309 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) -
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90
r4161 r5619 25 25 26 26 PUBLIC glob_sum ! used in many places 27 PUBLIC glob_sum_full ! used in many places 27 28 PUBLIC DDPDD ! also used in closea module 28 29 PUBLIC glob_min, glob_max … … 34 35 MODULE PROCEDURE glob_sum_1d, glob_sum_2d, glob_sum_3d, & 35 36 & glob_sum_2d_a, glob_sum_3d_a 37 END INTERFACE 38 INTERFACE glob_sum_full 39 MODULE PROCEDURE glob_sum_full_2d, glob_sum_full_3d 36 40 END INTERFACE 37 41 INTERFACE glob_min … … 156 160 ! 157 161 END FUNCTION glob_sum_3d_a 162 163 FUNCTION glob_sum_full_2d( ptab ) 164 !!---------------------------------------------------------------------- 165 !! *** FUNCTION glob_sum_2d *** 166 !! 167 !! ** Purpose : perform a sum in calling DDPDD routine 168 !!---------------------------------------------------------------------- 169 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptab 170 REAL(wp) :: glob_sum_full_2d ! global masked sum 171 !! 172 !!----------------------------------------------------------------------- 173 ! 174 glob_sum_full_2d = SUM( ptab(:,:) ) 175 IF( lk_mpp ) CALL mpp_sum( glob_sum_full_2d ) 176 ! 177 END FUNCTION glob_sum_full_2d 178 179 FUNCTION glob_sum_full_3d( ptab ) 180 !!---------------------------------------------------------------------- 181 !! *** FUNCTION glob_sum_3d *** 182 !! 183 !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine 184 !!---------------------------------------------------------------------- 185 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: ptab 186 REAL(wp) :: glob_sum_full_3d ! global masked sum 187 !! 188 INTEGER :: ji, jj, jk ! dummy loop indices 189 INTEGER :: ijpk ! local variables: size of ptab 190 !!----------------------------------------------------------------------- 191 ! 192 ijpk = SIZE(ptab,3) 193 ! 194 glob_sum_3d = 0.e0 195 DO jk = 1, ijpk 196 glob_sum_3d = glob_sum_3d + SUM( ptab(:,:,jk) ) 197 END DO 198 IF( lk_mpp ) CALL mpp_sum( glob_sum_full_3d ) 199 ! 200 END FUNCTION glob_sum_full_3d 201 158 202 159 203 #else … … 314 358 END FUNCTION glob_sum_3d_a 315 359 360 FUNCTION glob_sum_full_2d( ptab ) 361 !!---------------------------------------------------------------------- 362 !! *** FUNCTION glob_sum_2d *** 363 !! 364 !! ** Purpose : perform a sum in calling DDPDD routine 365 !!---------------------------------------------------------------------- 366 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptab 367 REAL(wp) :: glob_sum_full_2d ! global masked sum 368 !! 369 COMPLEX(wp):: ctmp 370 REAL(wp) :: ztmp 371 INTEGER :: ji, jj ! dummy loop indices 372 !!----------------------------------------------------------------------- 373 ! 374 ztmp = 0.e0 375 ctmp = CMPLX( 0.e0, 0.e0, wp ) 376 DO jj = 1, jpj 377 DO ji =1, jpi 378 ztmp = ptab(ji,jj) 379 CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 380 END DO 381 END DO 382 IF( lk_mpp ) CALL mpp_sum( ctmp ) ! sum over the global domain 383 glob_sum_full_2d = REAL(ctmp,wp) 384 ! 385 END FUNCTION glob_sum_full_2d 386 387 FUNCTION glob_sum_full_3d( ptab ) 388 !!---------------------------------------------------------------------- 389 !! *** FUNCTION glob_sum_3d *** 390 !! 391 !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine 392 !!---------------------------------------------------------------------- 393 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: ptab 394 REAL(wp) :: glob_sum_full_3d ! global masked sum 395 !! 396 COMPLEX(wp):: ctmp 397 REAL(wp) :: ztmp 398 INTEGER :: ji, jj, jk ! dummy loop indices 399 INTEGER :: ijpk ! local variables: size of ptab 400 !!----------------------------------------------------------------------- 401 ! 402 ijpk = SIZE(ptab,3) 403 ! 404 ztmp = 0.e0 405 ctmp = CMPLX( 0.e0, 0.e0, wp ) 406 DO jk = 1, ijpk 407 DO jj = 1, jpj 408 DO ji =1, jpi 409 ztmp = ptab(ji,jj,jk) 410 CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 411 END DO 412 END DO 413 END DO 414 IF( lk_mpp ) CALL mpp_sum( ctmp ) ! sum over the global domain 415 glob_sum_full_3d = REAL(ctmp,wp) 416 ! 417 END FUNCTION glob_sum_full_3d 418 419 420 316 421 #endif 317 422
Note: See TracChangeset
for help on using the changeset viewer.