Changeset 5260 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
- Timestamp:
- 2015-05-12T12:37:15+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
r4624 r5260 21 21 !! tra_dmp : update the tracer trend with the internal damping 22 22 !! tra_dmp_init : initialization, namlist read, parameters control 23 !! dtacof_zoom : restoring coefficient for zoom domain24 !! dtacof : restoring coefficient for global domain25 !! cofdis : compute the distance to the coastline26 23 !!---------------------------------------------------------------------- 27 24 USE oce ! ocean: variables 28 25 USE dom_oce ! ocean: domain variables 29 26 USE c1d ! 1D vertical configuration 30 USE trd mod_oce ! ocean: trendvariables31 USE trdtra ! active tracers: trends27 USE trd_oce ! trends: ocean variables 28 USE trdtra ! trends manager: tracers 32 29 USE zdf_oce ! ocean: vertical physics 33 30 USE phycst ! physical constants … … 39 36 USE wrk_nemo ! Memory allocation 40 37 USE timing ! Timing 38 USE iom 41 39 42 40 IMPLICIT NONE … … 45 43 PUBLIC tra_dmp ! routine called by step.F90 46 44 PUBLIC tra_dmp_init ! routine called by opa.F90 47 PUBLIC dtacof ! routine called by tradmp.F90, trcdmp.F90 and dyndmp.F9048 PUBLIC dtacof_zoom ! routine called by tradmp.F90, trcdmp.F90 and dyndmp.F9049 45 50 46 ! !!* Namelist namtra_dmp : T & S newtonian damping * 47 ! nn_zdmp and cn_resto are public as they are used by C1D/dyndmp.F90 51 48 LOGICAL , PUBLIC :: ln_tradmp !: internal damping flag 52 INTEGER , PUBLIC :: nn_hdmp ! = 0/-1/'latitude' for damping over T and S53 49 INTEGER , PUBLIC :: nn_zdmp ! = 0/1/2 flag for damping in the mixed layer 54 REAL(wp), PUBLIC :: rn_surf ! surface time scale for internal damping [days] 55 REAL(wp), PUBLIC :: rn_bot ! bottom time scale for internal damping [days] 56 REAL(wp), PUBLIC :: rn_dep ! depth of transition between rn_surf and rn_bot [meters] 57 INTEGER , PUBLIC :: nn_file ! = 1 create a damping.coeff NetCDF file 50 CHARACTER(LEN=200) , PUBLIC :: cn_resto ! name of netcdf file containing restoration coefficient field 51 ! 52 58 53 59 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: strdmp !: damping salinity trend (psu/s) … … 112 107 ! 113 108 CALL wrk_alloc( jpi, jpj, jpk, jpts, zts_dta ) 109 ! 114 110 ! !== input T-S data at kt ==! 115 111 CALL dta_tsd( kt, zts_dta ) ! read and interpolates T-S data at kt … … 172 168 ! 173 169 IF( l_trdtra ) THEN ! trend diagnostic 174 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_dmp, ttrdmp )175 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_dmp, strdmp )170 CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ttrdmp ) 171 CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, strdmp ) 176 172 ENDIF 177 173 ! ! Control print … … 194 190 !! ** Method : read the namtra_dmp namelist and check the parameters 195 191 !!---------------------------------------------------------------------- 196 NAMELIST/namtra_dmp/ ln_tradmp, nn_hdmp, nn_zdmp, rn_surf, rn_bot, rn_dep, nn_file 197 INTEGER :: ios ! Local integer output status for namelist read 198 !!---------------------------------------------------------------------- 199 200 REWIND( numnam_ref ) ! Namelist namtra_dmp in reference namelist : Temperature and salinity damping term 192 NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto 193 INTEGER :: ios ! Local integer for output status of namelist read 194 INTEGER :: imask ! File handle 195 !! 196 !!---------------------------------------------------------------------- 197 ! 198 REWIND( numnam_ref ) ! Namelist namtra_dmp in reference namelist : T & S relaxation 201 199 READ ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901) 202 200 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp ) 203 204 REWIND( numnam_cfg ) ! Namelist namtra_dmp in configuration namelist : Temperature and salinity damping term201 ! 202 REWIND( numnam_cfg ) ! Namelist namtra_dmp in configuration namelist : T & S relaxation 205 203 READ ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 ) 206 204 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp ) 207 205 IF(lwm) WRITE ( numond, namtra_dmp ) 208 209 IF( lzoom .AND. .NOT. lk_c1d ) nn_zdmp = 0 ! restoring to climatology at closed north or south boundaries 210 211 IF(lwp) THEN ! Namelist print 206 207 IF(lwp) THEN !Namelist print 212 208 WRITE(numout,*) 213 WRITE(numout,*) 'tra_dmp_init : T and S newtonian damping'209 WRITE(numout,*) 'tra_dmp_init : T and S newtonian relaxation' 214 210 WRITE(numout,*) '~~~~~~~' 215 WRITE(numout,*) ' Namelist namtra_dmp : set damping parameter' 216 WRITE(numout,*) ' add a damping term or not ln_tradmp = ', ln_tradmp 217 WRITE(numout,*) ' T and S damping option nn_hdmp = ', nn_hdmp 218 WRITE(numout,*) ' mixed layer damping option nn_zdmp = ', nn_zdmp, '(non-C1D zoom: forced to 0)' 219 WRITE(numout,*) ' surface time scale (days) rn_surf = ', rn_surf 220 WRITE(numout,*) ' bottom time scale (days) rn_bot = ', rn_bot 221 WRITE(numout,*) ' depth of transition (meters) rn_dep = ', rn_dep 222 WRITE(numout,*) ' create a damping.coeff file nn_file = ', nn_file 211 WRITE(numout,*) ' Namelist namtra_dmp : set relaxation parameters' 212 WRITE(numout,*) ' Apply relaxation or not ln_tradmp = ', ln_tradmp 213 WRITE(numout,*) ' mixed layer damping option nn_zdmp = ', nn_zdmp 214 WRITE(numout,*) ' Damping file name cn_resto = ', cn_resto 223 215 WRITE(numout,*) 224 216 ENDIF 225 217 226 IF( ln_tradmp ) THEN ! initialization for T-S damping 227 ! 218 IF( ln_tradmp) THEN 219 ! 220 !Allocate arrays 228 221 IF( tra_dmp_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' ) 229 ! 230 #if ! defined key_c1d 231 SELECT CASE ( nn_hdmp ) 232 CASE ( -1 ) ; IF(lwp) WRITE(numout,*) ' tracer damping in the Med & Red seas only' 233 CASE ( 1:90 ) ; IF(lwp) WRITE(numout,*) ' tracer damping poleward of', nn_hdmp, ' degrees' 234 CASE DEFAULT 235 WRITE(ctmp1,*) ' bad flag value for nn_hdmp = ', nn_hdmp 236 CALL ctl_stop(ctmp1) 222 223 !Check values of nn_zdmp 224 SELECT CASE (nn_zdmp) 225 CASE ( 0 ) ; IF(lwp) WRITE(numout,*) ' tracer damping as specified by mask' 226 CASE ( 1 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the turbocline' 227 CASE ( 2 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the mixed layer' 237 228 END SELECT 238 ! 239 #endif 240 SELECT CASE ( nn_zdmp ) 241 CASE ( 0 ) ; IF(lwp) WRITE(numout,*) ' tracer damping throughout the water column' 242 CASE ( 1 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the turbocline (avt > 5 cm2/s)' 243 CASE ( 2 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the mixed layer' 244 CASE DEFAULT 245 WRITE(ctmp1,*) 'bad flag value for nn_zdmp = ', nn_zdmp 246 CALL ctl_stop(ctmp1) 247 END SELECT 248 ! 229 230 !TG: Initialisation of dtatsd - Would it be better to have dmpdta routine 231 !so can damp to something other than intitial conditions files? 249 232 IF( .NOT.ln_tsd_tradmp ) THEN 250 233 CALL ctl_warn( 'tra_dmp_init: read T-S data not initialized, we force ln_tsd_tradmp=T' ) 251 234 CALL dta_tsd_init( ld_tradmp=ln_tradmp ) ! forces the initialisation of T-S data 252 235 ENDIF 253 ! 254 strdmp(:,:,:) = 0._wp ! internal damping salinity trend (used in asmtrj) 236 237 !initialise arrays - Are these actually used anywhere else? 238 strdmp(:,:,:) = 0._wp 255 239 ttrdmp(:,:,:) = 0._wp 256 ! ! Damping coefficients initialization 257 IF( lzoom .AND. .NOT. lk_c1d ) THEN ; CALL dtacof_zoom( resto )258 ELSE ; CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'TRA', resto)259 ENDIF260 !261 ENDIF262 ! 240 241 !Read in mask from file 242 CALL iom_open ( cn_resto, imask) 243 CALL iom_get ( imask, jpdom_autoglo, 'resto', resto) 244 CALL iom_close( imask ) 245 ENDIF 246 263 247 END SUBROUTINE tra_dmp_init 264 248 265 266 SUBROUTINE dtacof_zoom( presto )267 !!----------------------------------------------------------------------268 !! *** ROUTINE dtacof_zoom ***269 !!270 !! ** Purpose : Compute the damping coefficient for zoom domain271 !!272 !! ** Method : - set along closed boundary due to zoom a damping over273 !! 6 points with a max time scale of 5 days.274 !! - ORCA arctic/antarctic zoom: set the damping along275 !! south/north boundary over a latitude strip.276 !!277 !! ** Action : - resto, the damping coeff. for T and S278 !!----------------------------------------------------------------------279 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: presto ! restoring coeff. (s-1)280 !281 INTEGER :: ji, jj, jk, jn ! dummy loop indices282 REAL(wp) :: zlat, zlat0, zlat1, zlat2, z1_5d ! local scalar283 REAL(wp), DIMENSION(6) :: zfact ! 1Dworkspace284 !!----------------------------------------------------------------------285 !286 IF( nn_timing == 1 ) CALL timing_start( 'dtacof_zoom')287 !288 289 zfact(1) = 1._wp290 zfact(2) = 1._wp291 zfact(3) = 11._wp / 12._wp292 zfact(4) = 8._wp / 12._wp293 zfact(5) = 4._wp / 12._wp294 zfact(6) = 1._wp / 12._wp295 zfact(:) = zfact(:) / ( 5._wp * rday ) ! 5 days max restoring time scale296 297 presto(:,:,:) = 0._wp298 299 ! damping along the forced closed boundary over 6 grid-points300 DO jn = 1, 6301 IF( lzoom_w ) presto( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : ) = zfact(jn) ! west closed302 IF( lzoom_s ) presto( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : ) = zfact(jn) ! south closed303 IF( lzoom_e ) presto( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : ) = zfact(jn) ! east closed304 IF( lzoom_n ) presto( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : ) = zfact(jn) ! north closed305 END DO306 307 ! ! ====================================================308 IF( cp_cfz == "arctic" .OR. cp_cfz == "antarctic" ) THEN ! ORCA configuration : arctic or antarctic zoom309 ! ! ====================================================310 IF(lwp) WRITE(numout,*)311 IF(lwp .AND. cp_cfz == "arctic" ) WRITE(numout,*) ' dtacof_zoom : ORCA Arctic zoom'312 IF(lwp .AND. cp_cfz == "antarctic" ) WRITE(numout,*) ' dtacof_zoom : ORCA Antarctic zoom'313 IF(lwp) WRITE(numout,*)314 !315 ! ! Initialization :316 presto(:,:,:) = 0._wp317 zlat0 = 10._wp ! zlat0 : latitude strip where resto decreases318 zlat1 = 30._wp ! zlat1 : resto = 1 before zlat1319 zlat2 = zlat1 + zlat0 ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2320 z1_5d = 1._wp / ( 5._wp * rday ) ! z1_5d : 1 / 5days321 322 DO jk = 2, jpkm1 ! Compute arrays resto ; value for internal damping : 5 days323 DO jj = 1, jpj324 DO ji = 1, jpi325 zlat = ABS( gphit(ji,jj) )326 IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN327 presto(ji,jj,jk) = 0.5_wp * z1_5d * ( 1._wp - COS( rpi*(zlat2-zlat)/zlat0 ) )328 ELSEIF( zlat < zlat1 ) THEN329 presto(ji,jj,jk) = z1_5d330 ENDIF331 END DO332 END DO333 END DO334 !335 ENDIF336 ! ! Mask resto array337 presto(:,:,:) = presto(:,:,:) * tmask(:,:,:)338 !339 IF( nn_timing == 1 ) CALL timing_stop( 'dtacof_zoom')340 !341 END SUBROUTINE dtacof_zoom342 343 344 SUBROUTINE dtacof( kn_hdmp, pn_surf, pn_bot, pn_dep, &345 & kn_file, cdtype , presto )346 !!----------------------------------------------------------------------347 !! *** ROUTINE dtacof ***348 !!349 !! ** Purpose : Compute the damping coefficient350 !!351 !! ** Method : Arrays defining the damping are computed for each grid352 !! point for temperature and salinity (resto)353 !! Damping depends on distance to coast, depth and latitude354 !!355 !! ** Action : - resto, the damping coeff. for T and S356 !!----------------------------------------------------------------------357 USE iom358 USE ioipsl359 !!360 INTEGER , INTENT(in ) :: kn_hdmp ! damping option361 REAL(wp) , INTENT(in ) :: pn_surf ! surface time scale (days)362 REAL(wp) , INTENT(in ) :: pn_bot ! bottom time scale (days)363 REAL(wp) , INTENT(in ) :: pn_dep ! depth of transition (meters)364 INTEGER , INTENT(in ) :: kn_file ! save the damping coef on a file or not365 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA, TRC or DYN (tracer/dynamics indicator)366 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: presto ! restoring coeff. (s-1)367 !368 INTEGER :: ji, jj, jk ! dummy loop indices369 INTEGER :: ii0, ii1, ij0, ij1 ! local integers370 INTEGER :: inum0, icot ! - -371 REAL(wp) :: zinfl, zlon ! local scalars372 REAL(wp) :: zlat, zlat0, zlat1, zlat2 ! - -373 REAL(wp) :: zsdmp, zbdmp ! - -374 CHARACTER(len=20) :: cfile375 REAL(wp), POINTER, DIMENSION(: ) :: zhfac376 REAL(wp), POINTER, DIMENSION(:,: ) :: zmrs377 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdct378 !!----------------------------------------------------------------------379 !380 IF( nn_timing == 1 ) CALL timing_start('dtacof')381 !382 CALL wrk_alloc( jpk, zhfac )383 CALL wrk_alloc( jpi, jpj, zmrs )384 CALL wrk_alloc( jpi, jpj, jpk, zdct )385 #if defined key_c1d386 ! ! ====================387 ! ! C1D configuration : local domain388 ! ! ====================389 !390 IF(lwp) WRITE(numout,*)391 IF(lwp) WRITE(numout,*) ' dtacof : C1D 3x3 local domain'392 IF(lwp) WRITE(numout,*) ' -----------------------------'393 !394 presto(:,:,:) = 0._wp395 !396 zsdmp = 1._wp / ( pn_surf * rday )397 zbdmp = 1._wp / ( pn_bot * rday )398 DO jk = 2, jpkm1399 DO jj = 1, jpj400 DO ji = 1, jpi401 ! ONLY vertical variation from zsdmp (sea surface) to zbdmp (bottom)402 presto(ji,jj,jk) = zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep)403 END DO404 END DO405 END DO406 !407 presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:)408 #else409 ! ! ====================410 ! ! ORCA configuration : global domain411 ! ! ====================412 !413 IF(lwp) WRITE(numout,*)414 IF(lwp) WRITE(numout,*) ' dtacof : Global domain of ORCA'415 IF(lwp) WRITE(numout,*) ' ------------------------------'416 !417 presto(:,:,:) = 0._wp418 !419 IF( kn_hdmp > 0 ) THEN ! Damping poleward of 'nn_hdmp' degrees !420 ! !-----------------------------------------!421 IF(lwp) WRITE(numout,*)422 IF(lwp) WRITE(numout,*) ' Damping poleward of ', kn_hdmp, ' deg.'423 !424 CALL iom_open ( 'dist.coast.nc', icot, ldstop = .FALSE. )425 !426 IF( icot > 0 ) THEN ! distance-to-coast read in file427 CALL iom_get ( icot, jpdom_data, 'Tcoast', zdct )428 CALL iom_close( icot )429 ELSE ! distance-to-coast computed and saved in file (output in zdct)430 CALL cofdis( zdct )431 ENDIF432 433 ! ! Compute arrays resto434 zinfl = 1000.e3_wp ! distance of influence for damping term435 zlat0 = 10._wp ! latitude strip where resto decreases436 zlat1 = REAL( kn_hdmp ) ! resto = 0 between -zlat1 and zlat1437 zlat2 = zlat1 + zlat0 ! resto increases from 0 to 1 between |zlat1| and |zlat2|438 439 DO jj = 1, jpj440 DO ji = 1, jpi441 zlat = ABS( gphit(ji,jj) )442 IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN443 presto(ji,jj,1) = 0.5_wp * ( 1._wp - COS( rpi*(zlat-zlat1)/zlat0 ) )444 ELSEIF ( zlat > zlat2 ) THEN445 presto(ji,jj,1) = 1._wp446 ENDIF447 END DO448 END DO449 450 IF ( kn_hdmp == 20 ) THEN ! North Indian ocean (20N/30N x 45E/100E) : resto=0451 DO jj = 1, jpj452 DO ji = 1, jpi453 zlat = gphit(ji,jj)454 zlon = MOD( glamt(ji,jj), 360._wp )455 IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. 45._wp < zlon .AND. zlon < 100._wp ) THEN456 presto(ji,jj,1) = 0._wp457 ENDIF458 END DO459 END DO460 ENDIF461 462 zsdmp = 1._wp / ( pn_surf * rday )463 zbdmp = 1._wp / ( pn_bot * rday )464 DO jk = 2, jpkm1465 DO jj = 1, jpj466 DO ji = 1, jpi467 zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) )468 ! ... Decrease the value in the vicinity of the coast469 presto(ji,jj,jk) = presto(ji,jj,1 ) * 0.5_wp * ( 1._wp - COS( rpi*zdct(ji,jj,jk)/zinfl) )470 ! ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom)471 presto(ji,jj,jk) = presto(ji,jj,jk) * ( zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep) )472 END DO473 END DO474 END DO475 !476 ENDIF477 478 ! ! =========================479 ! ! Med and Red Sea damping (ORCA configuration only)480 ! ! =========================481 IF( cp_cfg == "orca" .AND. ( kn_hdmp > 0 .OR. kn_hdmp == -1 ) ) THEN482 IF(lwp)WRITE(numout,*)483 IF(lwp)WRITE(numout,*) ' ORCA configuration: Damping in Med and Red Seas'484 !485 zmrs(:,:) = 0._wp486 !487 SELECT CASE ( jp_cfg )488 ! ! =======================489 CASE ( 4 ) ! ORCA_R4 configuration490 ! ! =======================491 ij0 = 50 ; ij1 = 56 ! Mediterranean Sea492 493 ii0 = 81 ; ii1 = 91 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp494 ij0 = 50 ; ij1 = 55495 ii0 = 75 ; ii1 = 80 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp496 ij0 = 52 ; ij1 = 53497 ii0 = 70 ; ii1 = 74 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp498 ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea499 DO jk = 1, 17500 zhfac (jk) = 0.5_wp * ( 1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp ) ) / rday501 END DO502 DO jk = 18, jpkm1503 zhfac (jk) = 1._wp / rday504 END DO505 ! ! =======================506 CASE ( 2 ) ! ORCA_R2 configuration507 ! ! =======================508 ij0 = 96 ; ij1 = 110 ! Mediterranean Sea509 ii0 = 157 ; ii1 = 181 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp510 ij0 = 100 ; ij1 = 110511 ii0 = 144 ; ii1 = 156 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp512 ij0 = 100 ; ij1 = 103513 ii0 = 139 ; ii1 = 143 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp514 !515 ij0 = 101 ; ij1 = 102 ! Decrease before Gibraltar Strait516 ii0 = 139 ; ii1 = 141 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp517 ii0 = 142 ; ii1 = 142 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp518 ii0 = 143 ; ii1 = 143 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp519 ii0 = 144 ; ii1 = 144 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp520 !521 ij0 = 87 ; ij1 = 96 ! Red Sea522 ii0 = 147 ; ii1 = 163 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp523 !524 ij0 = 91 ; ij1 = 91 ! Decrease before Bab el Mandeb Strait525 ii0 = 153 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80_wp526 ij0 = 90 ; ij1 = 90527 ii0 = 153 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp528 ij0 = 89 ; ij1 = 89529 ii0 = 158 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp530 ij0 = 88 ; ij1 = 88531 ii0 = 160 ; ii1 = 163 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp532 ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea533 DO jk = 1, 17534 zhfac (jk) = 0.5_wp * ( 1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp ) ) / rday535 END DO536 DO jk = 18, jpkm1537 zhfac (jk) = 1._wp / rday538 END DO539 ! ! =======================540 CASE ( 05 ) ! ORCA_R05 configuration541 ! ! =======================542 ii0 = 568 ; ii1 = 574 ! Mediterranean Sea543 ij0 = 324 ; ij1 = 333 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp544 ii0 = 575 ; ii1 = 658545 ij0 = 314 ; ij1 = 366 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp546 !547 ii0 = 641 ; ii1 = 651 ! Black Sea (remaining part548 ij0 = 367 ; ij1 = 372 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp549 !550 ij0 = 324 ; ij1 = 333 ! Decrease before Gibraltar Strait551 ii0 = 565 ; ii1 = 565 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp552 ii0 = 566 ; ii1 = 566 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp553 ii0 = 567 ; ii1 = 567 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp554 !555 ii0 = 641 ; ii1 = 665 ! Red Sea556 ij0 = 270 ; ij1 = 310 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp557 !558 ii0 = 666 ; ii1 = 675 ! Decrease before Bab el Mandeb Strait559 ij0 = 270 ; ij1 = 290560 DO ji = mi0(ii0), mi1(ii1)561 zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1_wp * ABS( FLOAT(ji - mi1(ii1)) )562 END DO563 zsdmp = 1._wp / ( pn_surf * rday )564 zbdmp = 1._wp / ( pn_bot * rday )565 DO jk = 1, jpk566 zhfac(jk) = ( zbdmp + (zsdmp-zbdmp) * EXP( -fsdept(1,1,jk)/pn_dep ) )567 END DO568 ! ! ========================569 CASE ( 025 ) ! ORCA_R025 configuration570 ! ! ========================571 CALL ctl_stop( ' Not yet implemented in ORCA_R025' )572 !573 END SELECT574 575 DO jk = 1, jpkm1576 presto(:,:,jk) = zmrs(:,:) * zhfac(jk) + ( 1._wp - zmrs(:,:) ) * presto(:,:,jk)577 END DO578 579 ! Mask resto array and set to 0 first and last levels580 presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:)581 presto(:,:, 1 ) = 0._wp582 presto(:,:,jpk) = 0._wp583 ! !--------------------!584 ELSE ! No damping !585 ! !--------------------!586 CALL ctl_stop( 'Choose a correct value of nn_hdmp or put ln_tradmp to FALSE' )587 ENDIF588 #endif589 590 ! !--------------------------------!591 IF( kn_file == 1 ) THEN ! save damping coef. in a file !592 ! !--------------------------------!593 IF(lwp) WRITE(numout,*) ' create damping.coeff.nc file'594 IF( cdtype == 'TRA' ) cfile = 'damping.coeff'595 IF( cdtype == 'TRC' ) cfile = 'damping.coeff.trc'596 IF( cdtype == 'DYN' ) cfile = 'damping.coeff.dyn'597 cfile = TRIM( cfile )598 CALL iom_open ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )599 CALL iom_rstput( 0, 0, inum0, 'Resto', presto )600 CALL iom_close ( inum0 )601 ENDIF602 !603 CALL wrk_dealloc( jpk, zhfac)604 CALL wrk_dealloc( jpi, jpj, zmrs )605 CALL wrk_dealloc( jpi, jpj, jpk, zdct )606 !607 IF( nn_timing == 1 ) CALL timing_stop('dtacof')608 !609 END SUBROUTINE dtacof610 611 612 SUBROUTINE cofdis( pdct )613 !!----------------------------------------------------------------------614 !! *** ROUTINE cofdis ***615 !!616 !! ** Purpose : Compute the distance between ocean T-points and the617 !! ocean model coastlines. Save the distance in a NetCDF file.618 !!619 !! ** Method : For each model level, the distance-to-coast is620 !! computed as follows :621 !! - The coastline is defined as the serie of U-,V-,F-points622 !! that are at the ocean-land bound.623 !! - For each ocean T-point, the distance-to-coast is then624 !! computed as the smallest distance (on the sphere) between the625 !! T-point and all the coastline points.626 !! - For land T-points, the distance-to-coast is set to zero.627 !! C A U T I O N : Computation not yet implemented in mpp case.628 !!629 !! ** Action : - pdct, distance to the coastline (argument)630 !! - NetCDF file 'dist.coast.nc'631 !!----------------------------------------------------------------------632 USE ioipsl ! IOipsl librairy633 !!634 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: pdct ! distance to the coastline635 !!636 INTEGER :: ji, jj, jk, jl ! dummy loop indices637 INTEGER :: iju, ijt, icoast, itime, ierr, icot ! local integers638 CHARACTER (len=32) :: clname ! local name639 REAL(wp) :: zdate0 ! local scalar640 REAL(wp), POINTER, DIMENSION(:,:) :: zxt, zyt, zzt, zmask641 REAL(wp), POINTER, DIMENSION(: ) :: zxc, zyc, zzc, zdis ! temporary workspace642 LOGICAL , ALLOCATABLE, DIMENSION(:,:) :: llcotu, llcotv, llcotf ! 2D logical workspace643 !!----------------------------------------------------------------------644 !645 IF( nn_timing == 1 ) CALL timing_start('cofdis')646 !647 CALL wrk_alloc( jpi, jpj , zxt, zyt, zzt, zmask )648 CALL wrk_alloc( 3*jpi*jpj, zxc, zyc, zzc, zdis )649 ALLOCATE( llcotu(jpi,jpj), llcotv(jpi,jpj), llcotf(jpi,jpj) )650 !651 IF( lk_mpp ) CALL mpp_sum( ierr )652 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'cofdis: requested local arrays unavailable')653 654 ! 0. Initialization655 ! -----------------656 IF(lwp) WRITE(numout,*)657 IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline'658 IF(lwp) WRITE(numout,*) '~~~~~~'659 IF(lwp) WRITE(numout,*)660 IF( lk_mpp ) &661 & CALL ctl_stop(' Computation not yet implemented with key_mpp_...', &662 & ' Rerun the code on another computer or ', &663 & ' create the "dist.coast.nc" file using IDL' )664 665 pdct(:,:,:) = 0._wp666 zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) )667 zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) )668 zzt(:,:) = SIN( rad * gphit(:,:) )669 670 671 ! 1. Loop on vertical levels672 ! --------------------------673 ! ! ===============674 DO jk = 1, jpkm1 ! Horizontal slab675 ! ! ===============676 ! Define the coastline points (U, V and F)677 DO jj = 2, jpjm1678 DO ji = 2, jpim1679 zmask(ji,jj) = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &680 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )681 llcotu(ji,jj) = ( tmask(ji,jj, jk) + tmask(ji+1,jj ,jk) == 1._wp )682 llcotv(ji,jj) = ( tmask(ji,jj ,jk) + tmask(ji ,jj+1,jk) == 1._wp )683 llcotf(ji,jj) = ( zmask(ji,jj) > 0._wp ) .AND. ( zmask(ji,jj) < 4._wp )684 END DO685 END DO686 687 ! Lateral boundaries conditions688 llcotu(:, 1 ) = umask(:, 2 ,jk) == 1689 llcotu(:,jpj) = umask(:,jpjm1,jk) == 1690 llcotv(:, 1 ) = vmask(:, 2 ,jk) == 1691 llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1692 llcotf(:, 1 ) = fmask(:, 2 ,jk) == 1693 llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1694 695 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN696 llcotu( 1 ,:) = llcotu(jpim1,:)697 llcotu(jpi,:) = llcotu( 2 ,:)698 llcotv( 1 ,:) = llcotv(jpim1,:)699 llcotv(jpi,:) = llcotv( 2 ,:)700 llcotf( 1 ,:) = llcotf(jpim1,:)701 llcotf(jpi,:) = llcotf( 2 ,:)702 ELSE703 llcotu( 1 ,:) = umask( 2 ,:,jk) == 1704 llcotu(jpi,:) = umask(jpim1,:,jk) == 1705 llcotv( 1 ,:) = vmask( 2 ,:,jk) == 1706 llcotv(jpi,:) = vmask(jpim1,:,jk) == 1707 llcotf( 1 ,:) = fmask( 2 ,:,jk) == 1708 llcotf(jpi,:) = fmask(jpim1,:,jk) == 1709 ENDIF710 IF( nperio == 3 .OR. nperio == 4 ) THEN711 DO ji = 1, jpim1712 iju = jpi - ji + 1713 llcotu(ji,jpj ) = llcotu(iju,jpj-2)714 llcotf(ji,jpjm1) = llcotf(iju,jpj-2)715 llcotf(ji,jpj ) = llcotf(iju,jpj-3)716 END DO717 DO ji = jpi/2, jpim1718 iju = jpi - ji + 1719 llcotu(ji,jpjm1) = llcotu(iju,jpjm1)720 END DO721 DO ji = 2, jpi722 ijt = jpi - ji + 2723 llcotv(ji,jpjm1) = llcotv(ijt,jpj-2)724 llcotv(ji,jpj ) = llcotv(ijt,jpj-3)725 END DO726 ENDIF727 IF( nperio == 5 .OR. nperio == 6 ) THEN728 DO ji = 1, jpim1729 iju = jpi - ji730 llcotu(ji,jpj ) = llcotu(iju,jpjm1)731 llcotf(ji,jpj ) = llcotf(iju,jpj-2)732 END DO733 DO ji = jpi/2, jpim1734 iju = jpi - ji735 llcotf(ji,jpjm1) = llcotf(iju,jpjm1)736 END DO737 DO ji = 1, jpi738 ijt = jpi - ji + 1739 llcotv(ji,jpj ) = llcotv(ijt,jpjm1)740 END DO741 DO ji = jpi/2+1, jpi742 ijt = jpi - ji + 1743 llcotv(ji,jpjm1) = llcotv(ijt,jpjm1)744 END DO745 ENDIF746 747 ! Compute cartesian coordinates of coastline points748 ! and the number of coastline points749 icoast = 0750 DO jj = 1, jpj751 DO ji = 1, jpi752 IF( llcotf(ji,jj) ) THEN753 icoast = icoast + 1754 zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) )755 zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) )756 zzc(icoast) = SIN( rad*gphif(ji,jj) )757 ENDIF758 IF( llcotu(ji,jj) ) THEN759 icoast = icoast+1760 zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) )761 zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) )762 zzc(icoast) = SIN( rad*gphiu(ji,jj) )763 ENDIF764 IF( llcotv(ji,jj) ) THEN765 icoast = icoast+1766 zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) )767 zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) )768 zzc(icoast) = SIN( rad*gphiv(ji,jj) )769 ENDIF770 END DO771 END DO772 773 ! Distance for the T-points774 DO jj = 1, jpj775 DO ji = 1, jpi776 IF( tmask(ji,jj,jk) == 0._wp ) THEN777 pdct(ji,jj,jk) = 0._wp778 ELSE779 DO jl = 1, icoast780 zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2 &781 & + ( zyt(ji,jj) - zyc(jl) )**2 &782 & + ( zzt(ji,jj) - zzc(jl) )**2783 END DO784 pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) )785 ENDIF786 END DO787 END DO788 ! ! ===============789 END DO ! End of slab790 ! ! ===============791 792 793 ! 2. Create the distance to the coast file in NetCDF format794 ! ----------------------------------------------------------795 clname = 'dist.coast'796 itime = 0797 CALL ymds2ju( 0 , 1 , 1 , 0._wp , zdate0 )798 CALL restini( 'NONE', jpi , jpj , glamt, gphit , &799 & jpk , gdept_1d, clname, itime, zdate0, &800 & rdt , icot )801 CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct )802 CALL restclo( icot )803 !804 CALL wrk_dealloc( jpi, jpj , zxt, zyt, zzt, zmask )805 CALL wrk_dealloc( 3*jpi*jpj, zxc, zyc, zzc, zdis )806 DEALLOCATE( llcotu, llcotv, llcotf )807 !808 IF( nn_timing == 1 ) CALL timing_stop('cofdis')809 !810 END SUBROUTINE cofdis811 !!======================================================================812 249 END MODULE tradmp
Note: See TracChangeset
for help on using the changeset viewer.
