Changeset 7753 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
- Timestamp:
- 2017-03-03T12:46:59+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r7698 r7753 115 115 REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear 116 116 REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges 117 REAL(wp), POINTER, DIMENSION(:,:) :: z_ai118 117 ! 119 118 INTEGER, PARAMETER :: nitermax = 20 … … 123 122 IF( nn_timing == 1 ) CALL timing_start('limitd_me') 124 123 125 CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross , z_ai)124 CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross ) 126 125 127 126 ! conservation test … … 136 135 ! 137 136 138 !$OMP PARALLEL DO schedule(static) private(jj,ji)139 137 DO jj = 1, jpj ! Initialize arrays. 140 138 DO ji = 1, jpi … … 194 192 ! closing rate to a gross closing rate. 195 193 ! NOTE: 0 < aksum <= 1 196 !$OMP PARALLEL 197 !$OMP DO schedule(static) private(jj,ji) 198 DO jj = 1, jpj 199 DO ji = 1, jpi 200 closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj) 201 END DO 202 END DO 194 closing_gross(:,:) = closing_net(:,:) / aksum(:,:) 203 195 204 196 ! correction to closing rate and opening if closing rate is excessive … … 206 198 ! Reduce the closing rate if more than 100% of the open water 207 199 ! would be removed. Reduce the opening rate proportionately. 208 !$OMP DO schedule(static) private(jj,ji,za,zfac)209 200 DO jj = 1, jpj 210 201 DO ji = 1, jpi … … 225 216 ! would be removed. Reduce the opening rate proportionately. 226 217 DO jl = 1, jpl 227 !$OMP DO schedule(static) private(jj,ji,za,zfac)228 218 DO jj = 1, jpj 229 219 DO ji = 1, jpi … … 236 226 END DO 237 227 END DO 238 !$OMP END PARALLEL239 228 240 229 ! 3.3 Redistribute area, volume, and energy. … … 247 236 !-----------------------------------------------------------------------------! 248 237 ! This is in general not equal to one because of divergence during transport 249 !$OMP PARALLEL 250 !$OMP DO schedule(static) private(jj,ji) 251 DO jj = 1, jpj 252 DO ji = 1, jpi 253 asum(ji,jj) = 0._wp 254 z_ai(ji,jj) = 0._wp 255 END DO 256 END DO 257 DO jl = 1, jpl 258 !$OMP DO schedule(static) private(jj,ji) 259 DO jj = 1, jpj 260 DO ji = 1, jpi 261 z_ai(ji,jj) = z_ai(ji,jj) + a_i(ji,jj,jl) 262 END DO 263 END DO 264 END DO 265 !$OMP DO schedule(static) private(jj,ji) 266 DO jj = 1, jpj 267 DO ji = 1, jpi 268 asum(ji,jj) = ato_i(ji,jj) + z_ai(ji,jj) 269 END DO 270 END DO 238 asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 271 239 272 240 ! 3.5 Do we keep on iterating ??? … … 276 244 277 245 iterate_ridging = 0 278 !$OMP DO schedule(static) private(jj,ji)279 246 DO jj = 1, jpj 280 247 DO ji = 1, jpi … … 291 258 END DO 292 259 END DO 293 !$OMP END PARALLEL294 260 295 261 IF( lk_mpp ) CALL mpp_max( iterate_ridging ) … … 323 289 IF( ln_ctl ) CALL lim_prt3D( 'limitd_me' ) 324 290 325 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross , z_ai)291 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross ) 326 292 ! 327 293 IF( nn_timing == 1 ) CALL timing_stop('limitd_me') … … 340 306 REAL(wp) :: Gstari, astari, hrmean, zdummy ! local scalar 341 307 REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n 342 REAL(wp), POINTER, DIMENSION(:,:) :: z_ai343 308 !------------------------------------------------------------------------------! 344 309 345 310 CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 346 CALL wrk_alloc( jpi,jpj,z_ai )347 311 348 312 Gstari = 1.0/rn_gstar 349 313 astari = 1.0/rn_astar 350 !$OMP PARALLEL 351 !$OMP DO schedule(static) private(jj,ji) 352 DO jj = 1, jpj 353 DO ji = 1, jpi 354 aksum(ji,jj) = 0.0 355 END DO 356 END DO 357 !$OMP END DO NOWAIT 358 DO jl = 1, jpl 359 !$OMP DO schedule(static) private(jj,ji) 360 DO jj = 1, jpj 361 DO ji = 1, jpi 362 athorn(ji,jj,jl) = 0.0 363 aridge(ji,jj,jl) = 0.0 364 araft (ji,jj,jl) = 0.0 365 END DO 366 END DO 367 END DO 368 !$OMP END PARALLEL 314 aksum(:,:) = 0.0 315 athorn(:,:,:) = 0.0 316 aridge(:,:,:) = 0.0 317 araft (:,:,:) = 0.0 369 318 370 319 ! Zero out categories with very small areas 371 320 CALL lim_var_zapsmall 372 321 373 !$OMP PARALLEL374 322 ! Ice thickness needed for rafting 375 323 DO jl = 1, jpl 376 !$OMP DO schedule(static) private(jj,ji,rswitch)377 324 DO jj = 1, jpj 378 325 DO ji = 1, jpi … … 389 336 ! Compute total area of ice plus open water. 390 337 ! This is in general not equal to one because of divergence during transport 391 392 !$OMP DO schedule(static) private(jj,ji) 393 DO jj = 1, jpj 394 DO ji = 1, jpi 395 asum(ji,jj) = 0._wp 396 z_ai(ji,jj) = 0._wp 397 END DO 398 END DO 399 DO jl = 1, jpl 400 !$OMP DO schedule(static) private(jj,ji) 401 DO jj = 1, jpj 402 DO ji = 1, jpi 403 z_ai(ji,jj) = z_ai(ji,jj) + a_i(ji,jj,jl) 404 END DO 405 END DO 406 END DO 407 !$OMP DO schedule(static) private(jj,ji) 408 DO jj = 1, jpj 409 DO ji = 1, jpi 410 asum(ji,jj) = ato_i(ji,jj) + z_ai(ji,jj) 411 END DO 412 END DO 338 asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 339 413 340 ! Compute cumulative thickness distribution function 414 341 ! Compute the cumulative thickness distribution function Gsum, 415 342 ! where Gsum(n) is the fractional area in categories 0 to n. 416 343 ! initial value (in h = 0) equals open water area 417 !$OMP DO schedule(static) private(jj,ji) 418 DO jj = 1, jpj 419 DO ji = 1, jpi 420 Gsum(ji,jj,-1) = 0._wp 421 Gsum(ji,jj,0 ) = ato_i(ji,jj) 422 END DO 423 END DO 344 Gsum(:,:,-1) = 0._wp 345 Gsum(:,:,0 ) = ato_i(:,:) 424 346 ! for each value of h, you have to add ice concentration then 425 347 DO jl = 1, jpl 426 !$OMP DO schedule(static) private(jj,ji) 427 DO jj = 1, jpj 428 DO ji = 1, jpi 429 Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 430 END DO 431 END DO 348 Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 432 349 END DO 433 350 434 351 ! Normalize the cumulative distribution to 1 435 352 DO jl = 0, jpl 436 !$OMP DO schedule(static) private(jj,ji) 437 DO jj = 1, jpj 438 DO ji = 1, jpi 439 Gsum(ji,jj,jl) = Gsum(ji,jj,jl) / asum(ji,jj) 440 END DO 441 END DO 353 Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:) 442 354 END DO 443 !$OMP END PARALLEL444 355 445 356 ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) … … 458 369 IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 459 370 DO jl = 0, jpl 460 !$OMP PARALLEL DO schedule(static) private(jj,ji)461 371 DO jj = 1, jpj 462 372 DO ji = 1, jpi … … 477 387 ! 478 388 zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array 479 !$OMP PARALLEL480 389 DO jl = -1, jpl 481 !$OMP DO schedule(static) private(jj,ji) 482 DO jj = 1, jpj 483 DO ji = 1, jpi 484 Gsum(ji,jj,jl) = EXP( -Gsum(ji,jj,jl) * astari ) * zdummy 485 END DO 486 END DO 390 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 487 391 END DO 488 392 DO jl = 0, jpl 489 !$OMP DO schedule(static) private(jj,ji) 490 DO jj = 1, jpj 491 DO ji = 1, jpi 492 athorn(ji,jj,jl) = Gsum(ji,jj,jl-1) - Gsum(ji,jj,jl) 493 END DO 494 END DO 495 END DO 496 !$OMP END PARALLEL 393 athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 394 END DO 497 395 ! 498 396 ENDIF … … 502 400 ! 503 401 DO jl = 1, jpl 504 !$OMP PARALLEL DO schedule(static) private(jj,ji,zdummy)505 402 DO jj = 1, jpj 506 403 DO ji = 1, jpi … … 515 412 ! 516 413 DO jl = 1, jpl 517 !$OMP PARALLEL DO schedule(static) private(jj,ji) 518 DO jj = 1, jpj 519 DO ji = 1, jpi 520 aridge(ji,jj,jl) = athorn(ji,jj,jl) 521 END DO 522 END DO 414 aridge(:,:,jl) = athorn(:,:,jl) 523 415 END DO 524 416 ! … … 526 418 ! 527 419 DO jl = 1, jpl 528 !$OMP PARALLEL DO schedule(static) private(jj,ji) 529 DO jj = 1, jpj 530 DO ji = 1, jpi 531 araft(ji,jj,jl) = athorn(ji,jj,jl) 532 END DO 533 END DO 420 araft(:,:,jl) = athorn(:,:,jl) 534 421 END DO 535 422 ! … … 562 449 !----------------------------------------------------------------- 563 450 564 !$OMP PARALLEL 565 !$OMP DO schedule(static) private(jj,ji) 566 DO jj = 1, jpj 567 DO ji = 1, jpi 568 aksum(ji,jj) = athorn(ji,jj,0) 569 END DO 570 END DO 451 aksum(:,:) = athorn(:,:,0) 571 452 ! Transfer function 572 453 DO jl = 1, jpl !all categories have a specific transfer function 573 !$OMP DO schedule(static) private(jj,ji,hrmean)574 454 DO jj = 1, jpj 575 455 DO ji = 1, jpi … … 596 476 END DO 597 477 END DO 598 !$OMP END PARALLEL599 478 ! 600 479 CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 601 CALL wrk_dealloc( jpi,jpj,z_ai )602 480 ! 603 481 END SUBROUTINE lim_itd_me_ridgeprep … … 661 539 ! 1) Compute change in open water area due to closing and opening. 662 540 !------------------------------------------------------------------------------- 663 !$OMP PARALLEL DO schedule(static) private(jj,ji)664 541 DO jj = 1, jpj 665 542 DO ji = 1, jpi … … 691 568 END DO 692 569 693 !$OMP PARALLEL694 !$OMP DO schedule(static) private(ij,jj,ji)695 570 DO ij = 1, icells 696 571 ji = indxi(ij) ; jj = indxj(ij) … … 785 660 !-------------------------------------------------------------------- 786 661 DO jk = 1, nlay_i 787 !$OMP DO schedule(static) private(ij,jj,ji)788 662 DO ij = 1, icells 789 663 ji = indxi(ij) ; jj = indxj(ij) … … 813 687 DO jl2 = 1, jpl 814 688 ! over categories to which ridged/rafted ice is transferred 815 !$OMP DO schedule(static) private(ij,jj,ji,hL,hR,farea)816 689 DO ij = 1, icells 817 690 ji = indxi(ij) ; jj = indxj(ij) … … 848 721 ! Transfer ice energy to category jl2 by ridging 849 722 DO jk = 1, nlay_i 850 !$OMP DO schedule(static) private(ij,jj,ji)851 723 DO ij = 1, icells 852 724 ji = indxi(ij) ; jj = indxj(ij) … … 856 728 ! 857 729 END DO ! jl2 858 !$OMP END PARALLEL859 730 860 731 END DO ! jl1 (deforming categories) 732 861 733 ! 862 734 CALL wrk_dealloc( jpij, indxi, indxj ) … … 897 769 ! 1) Initialize 898 770 !------------------------------------------------------------------------------! 899 !$OMP PARALLEL DO schedule(static) private(jj,ji) 900 DO jj = 1, jpj 901 DO ji = 1, jpi 902 strength(ji,jj) = 0._wp 903 END DO 904 END DO 771 strength(:,:) = 0._wp 905 772 906 773 !------------------------------------------------------------------------------! … … 914 781 IF( kstrngth == 1 ) THEN 915 782 z1_3 = 1._wp / 3._wp 916 !$OMP PARALLEL917 783 DO jl = 1, jpl 918 !$OMP DO schedule(static) private(jj,ji)919 784 DO jj= 1, jpj 920 785 DO ji = 1, jpi … … 945 810 END DO 946 811 947 !$OMP DO schedule(static) private(jj,ji) 948 DO jj= 1, jpj 949 DO ji = 1, jpi 950 strength(ji,jj) = rn_pe_rdg * Cp * strength(ji,jj) / aksum(ji,jj) * tmask(ji,jj,1) 951 END DO 952 END DO 953 !$OMP END PARALLEL 812 strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1) 954 813 ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 955 814 ksmooth = 1 … … 959 818 !------------------------------------------------------------------------------! 960 819 ELSE ! kstrngth ne 1: Hibler (1979) form 961 !$OMP PARALLEL DO schedule(static) private(jj,ji) 962 DO jj= 1, jpj 963 DO ji = 1, jpi 964 ! 965 strength(ji,jj) = rn_pstar * vt_i(ji,jj) * EXP( - rn_crhg * ( 1._wp - at_i(ji,jj) ) ) * tmask(ji,jj,1) 966 END DO 967 END DO 820 ! 821 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) * tmask(:,:,1) 968 822 ! 969 823 ksmooth = 1 … … 976 830 ! CAN BE REMOVED 977 831 IF( ln_icestr_bvf ) THEN 978 !$OMP PARALLEL DO schedule(static) private(jj,ji)979 832 DO jj = 1, jpj 980 833 DO ji = 1, jpi … … 993 846 IF ( ksmooth == 1 ) THEN 994 847 995 !$OMP PARALLEL996 !$OMP DO schedule(static) private(jj,ji)997 848 DO jj = 2, jpjm1 998 849 DO ji = 2, jpim1 … … 1008 859 END DO 1009 860 1010 !$OMP DO schedule(static) private(jj,ji)1011 861 DO jj = 2, jpjm1 1012 862 DO ji = 2, jpim1 … … 1014 864 END DO 1015 865 END DO 1016 !$OMP END PARALLEL1017 866 CALL lbc_lnk( strength, 'T', 1. ) 1018 867 … … 1025 874 1026 875 IF ( numit == nit000 + nn_fsbc - 1 ) THEN 1027 !$OMP PARALLEL DO schedule(static) private(jj,ji) 1028 DO jj = 1, jpj 1029 DO ji = 1, jpi 1030 zstrp1(ji,jj) = 0._wp 1031 zstrp2(ji,jj) = 0._wp 1032 END DO 1033 END DO 876 zstrp1(:,:) = 0._wp 877 zstrp2(:,:) = 0._wp 1034 878 ENDIF 1035 879 1036 !$OMP PARALLEL DO schedule(static) private(jj,ji,numts_rm,zp)1037 880 DO jj = 2, jpjm1 1038 881 DO ji = 2, jpim1
Note: See TracChangeset
for help on using the changeset viewer.