- Timestamp:
- 2015-07-16T13:55:15+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
r4990 r5602 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 … … 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 50 !!gm why all namelist variable public???? only ln_tradmp should be sufficient51 45 52 46 ! !!* Namelist namtra_dmp : T & S newtonian damping * 47 ! nn_zdmp and cn_resto are public as they are used by C1D/dyndmp.F90 53 48 LOGICAL , PUBLIC :: ln_tradmp !: internal damping flag 54 INTEGER , PUBLIC :: nn_hdmp ! = 0/-1/'latitude' for damping over T and S55 49 INTEGER , PUBLIC :: nn_zdmp ! = 0/1/2 flag for damping in the mixed layer 56 REAL(wp), PUBLIC :: rn_surf ! surface time scale for internal damping [days] 57 REAL(wp), PUBLIC :: rn_bot ! bottom time scale for internal damping [days] 58 REAL(wp), PUBLIC :: rn_dep ! depth of transition between rn_surf and rn_bot [meters] 59 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 60 53 61 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: strdmp !: damping salinity trend (psu/s) … … 197 190 !! ** Method : read the namtra_dmp namelist and check the parameters 198 191 !!---------------------------------------------------------------------- 199 INTEGER :: ios ! Local integer output status for namelist read 200 !! 201 NAMELIST/namtra_dmp/ ln_tradmp, nn_hdmp, nn_zdmp, rn_surf, rn_bot, rn_dep, nn_file 202 !!---------------------------------------------------------------------- 203 ! 204 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 205 199 READ ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901) 206 200 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp ) 207 201 ! 208 REWIND( numnam_cfg ) ! Namelist namtra_dmp in configuration namelist : Temperature and salinity damping term202 REWIND( numnam_cfg ) ! Namelist namtra_dmp in configuration namelist : T & S relaxation 209 203 READ ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 ) 210 204 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp ) 211 205 IF(lwm) WRITE ( numond, namtra_dmp ) 212 213 IF( lzoom .AND. .NOT. lk_c1d ) nn_zdmp = 0 ! restoring to climatology at closed north or south boundaries 214 215 IF(lwp) THEN ! Namelist print 206 207 IF(lwp) THEN !Namelist print 216 208 WRITE(numout,*) 217 WRITE(numout,*) 'tra_dmp_init : T and S newtonian damping'209 WRITE(numout,*) 'tra_dmp_init : T and S newtonian relaxation' 218 210 WRITE(numout,*) '~~~~~~~' 219 WRITE(numout,*) ' Namelist namtra_dmp : set damping parameter' 220 WRITE(numout,*) ' add a damping term or not ln_tradmp = ', ln_tradmp 221 WRITE(numout,*) ' T and S damping option nn_hdmp = ', nn_hdmp 222 WRITE(numout,*) ' mixed layer damping option nn_zdmp = ', nn_zdmp, '(non-C1D zoom: forced to 0)' 223 WRITE(numout,*) ' surface time scale (days) rn_surf = ', rn_surf 224 WRITE(numout,*) ' bottom time scale (days) rn_bot = ', rn_bot 225 WRITE(numout,*) ' depth of transition (meters) rn_dep = ', rn_dep 226 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 227 215 WRITE(numout,*) 228 216 ENDIF 229 217 230 IF( ln_tradmp ) THEN ! initialization for T-S damping 231 ! 218 IF( ln_tradmp) THEN 219 ! 220 !Allocate arrays 232 221 IF( tra_dmp_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' ) 233 ! 234 !!gm I don't understand the specificities of c1d case...... 235 !!gm to be check with the autor of these lines 236 237 #if ! defined key_c1d 238 SELECT CASE ( nn_hdmp ) 239 CASE ( -1 ) ; IF(lwp) WRITE(numout,*) ' tracer damping in the Med & Red seas only' 240 CASE ( 1:90 ) ; IF(lwp) WRITE(numout,*) ' tracer damping poleward of', nn_hdmp, ' degrees' 241 CASE DEFAULT 242 WRITE(ctmp1,*) ' bad flag value for nn_hdmp = ', nn_hdmp 243 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' 244 228 END SELECT 245 ! 246 #endif 247 SELECT CASE ( nn_zdmp ) 248 CASE ( 0 ) ; IF(lwp) WRITE(numout,*) ' tracer damping throughout the water column' 249 CASE ( 1 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the turbocline (avt > 5 cm2/s)' 250 CASE ( 2 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the mixed layer' 251 CASE DEFAULT 252 WRITE(ctmp1,*) 'bad flag value for nn_zdmp = ', nn_zdmp 253 CALL ctl_stop(ctmp1) 254 END SELECT 255 ! 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? 256 232 IF( .NOT.ln_tsd_tradmp ) THEN 257 233 CALL ctl_warn( 'tra_dmp_init: read T-S data not initialized, we force ln_tsd_tradmp=T' ) 258 234 CALL dta_tsd_init( ld_tradmp=ln_tradmp ) ! forces the initialisation of T-S data 259 235 ENDIF 260 ! 261 strdmp(:,:,:) = 0._wp ! internal damping salinity trend (used in asmtrj) 236 237 !initialise arrays - Are these actually used anywhere else? 238 strdmp(:,:,:) = 0._wp 262 239 ttrdmp(:,:,:) = 0._wp 263 ! ! Damping coefficients initialization 264 IF( lzoom .AND. .NOT. lk_c1d ) THEN ; CALL dtacof_zoom( resto )265 ELSE ; CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'TRA', resto)266 ENDIF267 !268 ENDIF269 ! 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 270 247 END SUBROUTINE tra_dmp_init 271 248 272 273 SUBROUTINE dtacof_zoom( presto )274 !!----------------------------------------------------------------------275 !! *** ROUTINE dtacof_zoom ***276 !!277 !! ** Purpose : Compute the damping coefficient for zoom domain278 !!279 !! ** Method : - set along closed boundary due to zoom a damping over280 !! 6 points with a max time scale of 5 days.281 !! - ORCA arctic/antarctic zoom: set the damping along282 !! south/north boundary over a latitude strip.283 !!284 !! ** Action : - resto, the damping coeff. for T and S285 !!----------------------------------------------------------------------286 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: presto ! restoring coeff. (s-1)287 !288 INTEGER :: ji, jj, jk, jn ! dummy loop indices289 REAL(wp) :: zlat, zlat0, zlat1, zlat2, z1_5d ! local scalar290 REAL(wp), DIMENSION(6) :: zfact ! 1Dworkspace291 !!----------------------------------------------------------------------292 !293 IF( nn_timing == 1 ) CALL timing_start( 'dtacof_zoom')294 !295 296 zfact(1) = 1._wp297 zfact(2) = 1._wp298 zfact(3) = 11._wp / 12._wp299 zfact(4) = 8._wp / 12._wp300 zfact(5) = 4._wp / 12._wp301 zfact(6) = 1._wp / 12._wp302 zfact(:) = zfact(:) / ( 5._wp * rday ) ! 5 days max restoring time scale303 304 presto(:,:,:) = 0._wp305 306 ! damping along the forced closed boundary over 6 grid-points307 DO jn = 1, 6308 IF( lzoom_w ) presto( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : ) = zfact(jn) ! west closed309 IF( lzoom_s ) presto( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : ) = zfact(jn) ! south closed310 IF( lzoom_e ) presto( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : ) = zfact(jn) ! east closed311 IF( lzoom_n ) presto( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : ) = zfact(jn) ! north closed312 END DO313 314 ! ! ====================================================315 IF( cp_cfz == "arctic" .OR. cp_cfz == "antarctic" ) THEN ! ORCA configuration : arctic or antarctic zoom316 ! ! ====================================================317 IF(lwp) WRITE(numout,*)318 IF(lwp .AND. cp_cfz == "arctic" ) WRITE(numout,*) ' dtacof_zoom : ORCA Arctic zoom'319 IF(lwp .AND. cp_cfz == "antarctic" ) WRITE(numout,*) ' dtacof_zoom : ORCA Antarctic zoom'320 IF(lwp) WRITE(numout,*)321 !322 ! ! Initialization :323 presto(:,:,:) = 0._wp324 zlat0 = 10._wp ! zlat0 : latitude strip where resto decreases325 zlat1 = 30._wp ! zlat1 : resto = 1 before zlat1326 zlat2 = zlat1 + zlat0 ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2327 z1_5d = 1._wp / ( 5._wp * rday ) ! z1_5d : 1 / 5days328 329 DO jk = 2, jpkm1 ! Compute arrays resto ; value for internal damping : 5 days330 DO jj = 1, jpj331 DO ji = 1, jpi332 zlat = ABS( gphit(ji,jj) )333 IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN334 presto(ji,jj,jk) = 0.5_wp * z1_5d * ( 1._wp - COS( rpi*(zlat2-zlat)/zlat0 ) )335 ELSEIF( zlat < zlat1 ) THEN336 presto(ji,jj,jk) = z1_5d337 ENDIF338 END DO339 END DO340 END DO341 !342 ENDIF343 ! ! Mask resto array344 presto(:,:,:) = presto(:,:,:) * tmask(:,:,:)345 !346 IF( nn_timing == 1 ) CALL timing_stop( 'dtacof_zoom')347 !348 END SUBROUTINE dtacof_zoom349 350 351 SUBROUTINE dtacof( kn_hdmp, pn_surf, pn_bot, pn_dep, &352 & kn_file, cdtype , presto )353 !!----------------------------------------------------------------------354 !! *** ROUTINE dtacof ***355 !!356 !! ** Purpose : Compute the damping coefficient357 !!358 !! ** Method : Arrays defining the damping are computed for each grid359 !! point for temperature and salinity (resto)360 !! Damping depends on distance to coast, depth and latitude361 !!362 !! ** Action : - resto, the damping coeff. for T and S363 !!----------------------------------------------------------------------364 USE iom365 USE ioipsl366 !!367 INTEGER , INTENT(in ) :: kn_hdmp ! damping option368 REAL(wp) , INTENT(in ) :: pn_surf ! surface time scale (days)369 REAL(wp) , INTENT(in ) :: pn_bot ! bottom time scale (days)370 REAL(wp) , INTENT(in ) :: pn_dep ! depth of transition (meters)371 INTEGER , INTENT(in ) :: kn_file ! save the damping coef on a file or not372 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA, TRC or DYN (tracer/dynamics indicator)373 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: presto ! restoring coeff. (s-1)374 !375 INTEGER :: ji, jj, jk ! dummy loop indices376 INTEGER :: ii0, ii1, ij0, ij1 ! local integers377 INTEGER :: inum0, icot ! - -378 REAL(wp) :: zinfl, zlon ! local scalars379 REAL(wp) :: zlat, zlat0, zlat1, zlat2 ! - -380 REAL(wp) :: zsdmp, zbdmp ! - -381 CHARACTER(len=20) :: cfile382 REAL(wp), POINTER, DIMENSION(: ) :: zhfac383 REAL(wp), POINTER, DIMENSION(:,: ) :: zmrs384 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdct385 !!----------------------------------------------------------------------386 !387 IF( nn_timing == 1 ) CALL timing_start('dtacof')388 !389 CALL wrk_alloc( jpk, zhfac )390 CALL wrk_alloc( jpi, jpj, zmrs )391 CALL wrk_alloc( jpi, jpj, jpk, zdct )392 #if defined key_c1d393 ! ! ====================394 ! ! C1D configuration : local domain395 ! ! ====================396 !397 IF(lwp) WRITE(numout,*)398 IF(lwp) WRITE(numout,*) ' dtacof : C1D 3x3 local domain'399 IF(lwp) WRITE(numout,*) ' -----------------------------'400 !401 presto(:,:,:) = 0._wp402 !403 zsdmp = 1._wp / ( pn_surf * rday )404 zbdmp = 1._wp / ( pn_bot * rday )405 DO jk = 2, jpkm1406 DO jj = 1, jpj407 DO ji = 1, jpi408 ! ONLY vertical variation from zsdmp (sea surface) to zbdmp (bottom)409 presto(ji,jj,jk) = zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep)410 END DO411 END DO412 END DO413 !414 presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:)415 #else416 ! ! ====================417 ! ! ORCA configuration : global domain418 ! ! ====================419 !420 IF(lwp) WRITE(numout,*)421 IF(lwp) WRITE(numout,*) ' dtacof : Global domain of ORCA'422 IF(lwp) WRITE(numout,*) ' ------------------------------'423 !424 presto(:,:,:) = 0._wp425 !426 IF( kn_hdmp > 0 ) THEN ! Damping poleward of 'nn_hdmp' degrees !427 ! !-----------------------------------------!428 IF(lwp) WRITE(numout,*)429 IF(lwp) WRITE(numout,*) ' Damping poleward of ', kn_hdmp, ' deg.'430 !431 CALL iom_open ( 'dist.coast.nc', icot, ldstop = .FALSE. )432 !433 IF( icot > 0 ) THEN ! distance-to-coast read in file434 CALL iom_get ( icot, jpdom_data, 'Tcoast', zdct )435 CALL iom_close( icot )436 ELSE ! distance-to-coast computed and saved in file (output in zdct)437 CALL cofdis( zdct )438 ENDIF439 440 ! ! Compute arrays resto441 zinfl = 1000.e3_wp ! distance of influence for damping term442 zlat0 = 10._wp ! latitude strip where resto decreases443 zlat1 = REAL( kn_hdmp ) ! resto = 0 between -zlat1 and zlat1444 zlat2 = zlat1 + zlat0 ! resto increases from 0 to 1 between |zlat1| and |zlat2|445 446 DO jj = 1, jpj447 DO ji = 1, jpi448 zlat = ABS( gphit(ji,jj) )449 IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN450 presto(ji,jj,1) = 0.5_wp * ( 1._wp - COS( rpi*(zlat-zlat1)/zlat0 ) )451 ELSEIF ( zlat > zlat2 ) THEN452 presto(ji,jj,1) = 1._wp453 ENDIF454 END DO455 END DO456 457 IF ( kn_hdmp == 20 ) THEN ! North Indian ocean (20N/30N x 45E/100E) : resto=0458 DO jj = 1, jpj459 DO ji = 1, jpi460 zlat = gphit(ji,jj)461 zlon = MOD( glamt(ji,jj), 360._wp )462 IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. 45._wp < zlon .AND. zlon < 100._wp ) THEN463 presto(ji,jj,1) = 0._wp464 ENDIF465 END DO466 END DO467 ENDIF468 469 zsdmp = 1._wp / ( pn_surf * rday )470 zbdmp = 1._wp / ( pn_bot * rday )471 DO jk = 2, jpkm1472 DO jj = 1, jpj473 DO ji = 1, jpi474 zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) )475 ! ... Decrease the value in the vicinity of the coast476 presto(ji,jj,jk) = presto(ji,jj,1 ) * 0.5_wp * ( 1._wp - COS( rpi*zdct(ji,jj,jk)/zinfl) )477 ! ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom)478 presto(ji,jj,jk) = presto(ji,jj,jk) * ( zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep) )479 END DO480 END DO481 END DO482 !483 ENDIF484 485 ! ! =========================486 ! ! Med and Red Sea damping (ORCA configuration only)487 ! ! =========================488 IF( cp_cfg == "orca" .AND. ( kn_hdmp > 0 .OR. kn_hdmp == -1 ) ) THEN489 IF(lwp)WRITE(numout,*)490 IF(lwp)WRITE(numout,*) ' ORCA configuration: Damping in Med and Red Seas'491 !492 zmrs(:,:) = 0._wp493 !494 SELECT CASE ( jp_cfg )495 ! ! =======================496 CASE ( 4 ) ! ORCA_R4 configuration497 ! ! =======================498 ij0 = 50 ; ij1 = 56 ! Mediterranean Sea499 500 ii0 = 81 ; ii1 = 91 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp501 ij0 = 50 ; ij1 = 55502 ii0 = 75 ; ii1 = 80 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp503 ij0 = 52 ; ij1 = 53504 ii0 = 70 ; ii1 = 74 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp505 ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea506 DO jk = 1, 17507 zhfac (jk) = 0.5_wp * ( 1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp ) ) / rday508 END DO509 DO jk = 18, jpkm1510 zhfac (jk) = 1._wp / rday511 END DO512 ! ! =======================513 CASE ( 2 ) ! ORCA_R2 configuration514 ! ! =======================515 ij0 = 96 ; ij1 = 110 ! Mediterranean Sea516 ii0 = 157 ; ii1 = 181 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp517 ij0 = 100 ; ij1 = 110518 ii0 = 144 ; ii1 = 156 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp519 ij0 = 100 ; ij1 = 103520 ii0 = 139 ; ii1 = 143 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp521 !522 ij0 = 101 ; ij1 = 102 ! Decrease before Gibraltar Strait523 ii0 = 139 ; ii1 = 141 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp524 ii0 = 142 ; ii1 = 142 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp525 ii0 = 143 ; ii1 = 143 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp526 ii0 = 144 ; ii1 = 144 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp527 !528 ij0 = 87 ; ij1 = 96 ! Red Sea529 ii0 = 147 ; ii1 = 163 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp530 !531 ij0 = 91 ; ij1 = 91 ! Decrease before Bab el Mandeb Strait532 ii0 = 153 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80_wp533 ij0 = 90 ; ij1 = 90534 ii0 = 153 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp535 ij0 = 89 ; ij1 = 89536 ii0 = 158 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp537 ij0 = 88 ; ij1 = 88538 ii0 = 160 ; ii1 = 163 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp539 ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea540 DO jk = 1, 17541 zhfac (jk) = 0.5_wp * ( 1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp ) ) / rday542 END DO543 DO jk = 18, jpkm1544 zhfac (jk) = 1._wp / rday545 END DO546 ! ! =======================547 CASE ( 05 ) ! ORCA_R05 configuration548 ! ! =======================549 ii0 = 568 ; ii1 = 574 ! Mediterranean Sea550 ij0 = 324 ; ij1 = 333 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp551 ii0 = 575 ; ii1 = 658552 ij0 = 314 ; ij1 = 366 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp553 !554 ii0 = 641 ; ii1 = 651 ! Black Sea (remaining part555 ij0 = 367 ; ij1 = 372 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp556 !557 ij0 = 324 ; ij1 = 333 ! Decrease before Gibraltar Strait558 ii0 = 565 ; ii1 = 565 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp559 ii0 = 566 ; ii1 = 566 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp560 ii0 = 567 ; ii1 = 567 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp561 !562 ii0 = 641 ; ii1 = 665 ! Red Sea563 ij0 = 270 ; ij1 = 310 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp564 !565 ii0 = 666 ; ii1 = 675 ! Decrease before Bab el Mandeb Strait566 ij0 = 270 ; ij1 = 290567 DO ji = mi0(ii0), mi1(ii1)568 zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1_wp * ABS( FLOAT(ji - mi1(ii1)) )569 END DO570 zsdmp = 1._wp / ( pn_surf * rday )571 zbdmp = 1._wp / ( pn_bot * rday )572 DO jk = 1, jpk573 zhfac(jk) = ( zbdmp + (zsdmp-zbdmp) * EXP( -fsdept(1,1,jk)/pn_dep ) )574 END DO575 ! ! ========================576 CASE ( 025 ) ! ORCA_R025 configuration577 ! ! ========================578 CALL ctl_stop( ' Not yet implemented in ORCA_R025' )579 !580 END SELECT581 582 DO jk = 1, jpkm1583 presto(:,:,jk) = zmrs(:,:) * zhfac(jk) + ( 1._wp - zmrs(:,:) ) * presto(:,:,jk)584 END DO585 586 ! Mask resto array and set to 0 first and last levels587 presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:)588 presto(:,:, 1 ) = 0._wp589 presto(:,:,jpk) = 0._wp590 ! !--------------------!591 ELSE ! No damping !592 ! !--------------------!593 CALL ctl_stop( 'Choose a correct value of nn_hdmp or put ln_tradmp to FALSE' )594 ENDIF595 #endif596 597 ! !--------------------------------!598 IF( kn_file == 1 ) THEN ! save damping coef. in a file !599 ! !--------------------------------!600 IF(lwp) WRITE(numout,*) ' create damping.coeff.nc file'601 IF( cdtype == 'TRA' ) cfile = 'damping.coeff'602 IF( cdtype == 'TRC' ) cfile = 'damping.coeff.trc'603 IF( cdtype == 'DYN' ) cfile = 'damping.coeff.dyn'604 cfile = TRIM( cfile )605 CALL iom_open ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )606 CALL iom_rstput( 0, 0, inum0, 'Resto', presto )607 CALL iom_close ( inum0 )608 ENDIF609 !610 CALL wrk_dealloc( jpk, zhfac)611 CALL wrk_dealloc( jpi, jpj, zmrs )612 CALL wrk_dealloc( jpi, jpj, jpk, zdct )613 !614 IF( nn_timing == 1 ) CALL timing_stop('dtacof')615 !616 END SUBROUTINE dtacof617 618 619 SUBROUTINE cofdis( pdct )620 !!----------------------------------------------------------------------621 !! *** ROUTINE cofdis ***622 !!623 !! ** Purpose : Compute the distance between ocean T-points and the624 !! ocean model coastlines. Save the distance in a NetCDF file.625 !!626 !! ** Method : For each model level, the distance-to-coast is627 !! computed as follows :628 !! - The coastline is defined as the serie of U-,V-,F-points629 !! that are at the ocean-land bound.630 !! - For each ocean T-point, the distance-to-coast is then631 !! computed as the smallest distance (on the sphere) between the632 !! T-point and all the coastline points.633 !! - For land T-points, the distance-to-coast is set to zero.634 !! C A U T I O N : Computation not yet implemented in mpp case.635 !!636 !! ** Action : - pdct, distance to the coastline (argument)637 !! - NetCDF file 'dist.coast.nc'638 !!----------------------------------------------------------------------639 USE ioipsl ! IOipsl librairy640 !!641 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: pdct ! distance to the coastline642 !!643 INTEGER :: ji, jj, jk, jl ! dummy loop indices644 INTEGER :: iju, ijt, icoast, itime, ierr, icot ! local integers645 CHARACTER (len=32) :: clname ! local name646 REAL(wp) :: zdate0 ! local scalar647 REAL(wp), POINTER, DIMENSION(:,:) :: zxt, zyt, zzt, zmask648 REAL(wp), POINTER, DIMENSION(: ) :: zxc, zyc, zzc, zdis ! temporary workspace649 LOGICAL , ALLOCATABLE, DIMENSION(:,:) :: llcotu, llcotv, llcotf ! 2D logical workspace650 !!----------------------------------------------------------------------651 !652 IF( nn_timing == 1 ) CALL timing_start('cofdis')653 !654 CALL wrk_alloc( jpi, jpj , zxt, zyt, zzt, zmask )655 CALL wrk_alloc( 3*jpi*jpj, zxc, zyc, zzc, zdis )656 ALLOCATE( llcotu(jpi,jpj), llcotv(jpi,jpj), llcotf(jpi,jpj) )657 !658 IF( lk_mpp ) CALL mpp_sum( ierr )659 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'cofdis: requested local arrays unavailable')660 661 ! 0. Initialization662 ! -----------------663 IF(lwp) WRITE(numout,*)664 IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline'665 IF(lwp) WRITE(numout,*) '~~~~~~'666 IF(lwp) WRITE(numout,*)667 IF( lk_mpp ) &668 & CALL ctl_stop(' Computation not yet implemented with key_mpp_...', &669 & ' Rerun the code on another computer or ', &670 & ' create the "dist.coast.nc" file using IDL' )671 672 pdct(:,:,:) = 0._wp673 zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) )674 zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) )675 zzt(:,:) = SIN( rad * gphit(:,:) )676 677 678 ! 1. Loop on vertical levels679 ! --------------------------680 ! ! ===============681 DO jk = 1, jpkm1 ! Horizontal slab682 ! ! ===============683 ! Define the coastline points (U, V and F)684 DO jj = 2, jpjm1685 DO ji = 2, jpim1686 zmask(ji,jj) = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &687 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )688 llcotu(ji,jj) = ( tmask(ji,jj, jk) + tmask(ji+1,jj ,jk) == 1._wp )689 llcotv(ji,jj) = ( tmask(ji,jj ,jk) + tmask(ji ,jj+1,jk) == 1._wp )690 llcotf(ji,jj) = ( zmask(ji,jj) > 0._wp ) .AND. ( zmask(ji,jj) < 4._wp )691 END DO692 END DO693 694 ! Lateral boundaries conditions695 llcotu(:, 1 ) = umask(:, 2 ,jk) == 1696 llcotu(:,jpj) = umask(:,jpjm1,jk) == 1697 llcotv(:, 1 ) = vmask(:, 2 ,jk) == 1698 llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1699 llcotf(:, 1 ) = fmask(:, 2 ,jk) == 1700 llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1701 702 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN703 llcotu( 1 ,:) = llcotu(jpim1,:)704 llcotu(jpi,:) = llcotu( 2 ,:)705 llcotv( 1 ,:) = llcotv(jpim1,:)706 llcotv(jpi,:) = llcotv( 2 ,:)707 llcotf( 1 ,:) = llcotf(jpim1,:)708 llcotf(jpi,:) = llcotf( 2 ,:)709 ELSE710 llcotu( 1 ,:) = umask( 2 ,:,jk) == 1711 llcotu(jpi,:) = umask(jpim1,:,jk) == 1712 llcotv( 1 ,:) = vmask( 2 ,:,jk) == 1713 llcotv(jpi,:) = vmask(jpim1,:,jk) == 1714 llcotf( 1 ,:) = fmask( 2 ,:,jk) == 1715 llcotf(jpi,:) = fmask(jpim1,:,jk) == 1716 ENDIF717 IF( nperio == 3 .OR. nperio == 4 ) THEN718 DO ji = 1, jpim1719 iju = jpi - ji + 1720 llcotu(ji,jpj ) = llcotu(iju,jpj-2)721 llcotf(ji,jpjm1) = llcotf(iju,jpj-2)722 llcotf(ji,jpj ) = llcotf(iju,jpj-3)723 END DO724 DO ji = jpi/2, jpim1725 iju = jpi - ji + 1726 llcotu(ji,jpjm1) = llcotu(iju,jpjm1)727 END DO728 DO ji = 2, jpi729 ijt = jpi - ji + 2730 llcotv(ji,jpjm1) = llcotv(ijt,jpj-2)731 llcotv(ji,jpj ) = llcotv(ijt,jpj-3)732 END DO733 ENDIF734 IF( nperio == 5 .OR. nperio == 6 ) THEN735 DO ji = 1, jpim1736 iju = jpi - ji737 llcotu(ji,jpj ) = llcotu(iju,jpjm1)738 llcotf(ji,jpj ) = llcotf(iju,jpj-2)739 END DO740 DO ji = jpi/2, jpim1741 iju = jpi - ji742 llcotf(ji,jpjm1) = llcotf(iju,jpjm1)743 END DO744 DO ji = 1, jpi745 ijt = jpi - ji + 1746 llcotv(ji,jpj ) = llcotv(ijt,jpjm1)747 END DO748 DO ji = jpi/2+1, jpi749 ijt = jpi - ji + 1750 llcotv(ji,jpjm1) = llcotv(ijt,jpjm1)751 END DO752 ENDIF753 754 ! Compute cartesian coordinates of coastline points755 ! and the number of coastline points756 icoast = 0757 DO jj = 1, jpj758 DO ji = 1, jpi759 IF( llcotf(ji,jj) ) THEN760 icoast = icoast + 1761 zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) )762 zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) )763 zzc(icoast) = SIN( rad*gphif(ji,jj) )764 ENDIF765 IF( llcotu(ji,jj) ) THEN766 icoast = icoast+1767 zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) )768 zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) )769 zzc(icoast) = SIN( rad*gphiu(ji,jj) )770 ENDIF771 IF( llcotv(ji,jj) ) THEN772 icoast = icoast+1773 zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) )774 zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) )775 zzc(icoast) = SIN( rad*gphiv(ji,jj) )776 ENDIF777 END DO778 END DO779 780 ! Distance for the T-points781 DO jj = 1, jpj782 DO ji = 1, jpi783 IF( tmask(ji,jj,jk) == 0._wp ) THEN784 pdct(ji,jj,jk) = 0._wp785 ELSE786 DO jl = 1, icoast787 zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2 &788 & + ( zyt(ji,jj) - zyc(jl) )**2 &789 & + ( zzt(ji,jj) - zzc(jl) )**2790 END DO791 pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) )792 ENDIF793 END DO794 END DO795 ! ! ===============796 END DO ! End of slab797 ! ! ===============798 799 800 ! 2. Create the distance to the coast file in NetCDF format801 ! ----------------------------------------------------------802 clname = 'dist.coast'803 itime = 0804 CALL ymds2ju( 0 , 1 , 1 , 0._wp , zdate0 )805 CALL restini( 'NONE', jpi , jpj , glamt, gphit , &806 & jpk , gdept_1d, clname, itime, zdate0, &807 & rdt , icot )808 CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct )809 CALL restclo( icot )810 !811 CALL wrk_dealloc( jpi, jpj , zxt, zyt, zzt, zmask )812 CALL wrk_dealloc( 3*jpi*jpj, zxc, zyc, zzc, zdis )813 DEALLOCATE( llcotu, llcotv, llcotf )814 !815 IF( nn_timing == 1 ) CALL timing_stop('cofdis')816 !817 END SUBROUTINE cofdis818 !!======================================================================819 249 END MODULE tradmp
Note: See TracChangeset
for help on using the changeset viewer.