Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/ASM/asminc.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/ASM/asminc.F90 (modified) (54 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/ASM/asminc.F90
r13295 r14789 9 9 !! ! 2007-04 (A. Weaver) Merge with OPAVAR/NEMOVAR 10 10 !! NEMO 3.3 ! 2010-05 (D. Lea) Update to work with NEMO v3.2 11 !! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init 11 !! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init 12 12 !! 3.4 ! 2012-10 (A. Weaver and K. Mogensen) Fix for direct initialization 13 13 !! ! 2014-09 (D. Lea) Local calc_date removed use routine from OBS … … 26 26 USE par_oce ! Ocean space and time domain variables 27 27 USE dom_oce ! Ocean space and time domain 28 USE domtile 28 29 USE domvvl ! domain: variable volume level 29 30 USE ldfdyn ! lateral diffusion: eddy viscosity coefficients … … 31 32 USE zpshde ! Partial step : Horizontal Derivative 32 33 USE asmpar ! Parameters for the assmilation interface 33 USE asmbkg ! 34 USE asmbkg ! 34 35 USE c1d ! 1D initialization 35 36 USE sbc_oce ! Surface boundary condition variables. … … 45 46 IMPLICIT NONE 46 47 PRIVATE 47 48 48 49 PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights 49 50 PUBLIC tra_asm_inc !: Apply the tracer (T and S) increments … … 72 73 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkg , v_bkg !: Background u- & v- velocity components 73 74 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkginc, s_bkginc !: Increment to the background T & S 74 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components 75 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components 75 76 REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step 76 77 #if defined key_asminc … … 80 81 INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term 81 82 INTEGER , PUBLIC :: nitdin !: Time step of the background state for direct initialization 82 INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval 83 INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval 83 84 INTEGER , PUBLIC :: nitiaufin !: Time step of the end of the IAU interval 84 ! 85 ! 85 86 INTEGER , PUBLIC :: niaufn !: Type of IAU weighing function: = 0 Constant weighting 86 ! !: = 1 Linear hat-like, centred in middle of IAU interval 87 ! !: = 1 Linear hat-like, centred in middle of IAU interval 87 88 REAL(wp), PUBLIC :: salfixmin !: Ensure that the salinity is larger than this value if (ln_salfix) 88 89 … … 106 107 !!---------------------------------------------------------------------- 107 108 !! *** ROUTINE asm_inc_init *** 108 !! 109 !! 109 110 !! ** Purpose : Initialize the assimilation increment and IAU weights. 110 111 !! 111 112 !! ** Method : Initialize the assimilation increment and IAU weights. 112 113 !! 113 !! ** Action : 114 !! ** Action : 114 115 !!---------------------------------------------------------------------- 115 116 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices … … 263 264 ! 264 265 ! !--------------------------------------------------------- 265 IF( niaufn == 0 ) THEN ! Constant IAU forcing 266 IF( niaufn == 0 ) THEN ! Constant IAU forcing 266 267 ! !--------------------------------------------------------- 267 268 DO jt = 1, iiauper … … 269 270 END DO 270 271 ! !--------------------------------------------------------- 271 ELSEIF ( niaufn == 1 ) THEN ! Linear hat-like, centred in middle of IAU interval 272 ELSEIF ( niaufn == 1 ) THEN ! Linear hat-like, centred in middle of IAU interval 272 273 ! !--------------------------------------------------------- 273 274 ! Compute the normalization factor 274 275 znorm = 0._wp 275 276 IF( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval 276 imid = iiauper / 2 277 imid = iiauper / 2 277 278 DO jt = 1, imid 278 279 znorm = znorm + REAL( jt ) … … 280 281 znorm = 2.0 * znorm 281 282 ELSE ! Odd number of time steps in IAU interval 282 imid = ( iiauper + 1 ) / 2 283 imid = ( iiauper + 1 ) / 2 283 284 DO jt = 1, imid - 1 284 285 znorm = znorm + REAL( jt ) … … 307 308 DO jt = 1, icycper 308 309 ztotwgt = ztotwgt + wgtiau(jt) 309 WRITE(numout,*) ' ', jt, ' ', wgtiau(jt) 310 END DO 310 WRITE(numout,*) ' ', jt, ' ', wgtiau(jt) 311 END DO 311 312 WRITE(numout,*) ' ===================================' 312 313 WRITE(numout,*) ' Time-integrated weight = ', ztotwgt 313 314 WRITE(numout,*) ' ===================================' 314 315 ENDIF 315 316 316 317 ENDIF 317 318 … … 338 339 CALL iom_open( c_asminc, inum ) 339 340 ! 340 CALL iom_get( inum, 'time' , zdate_inc ) 341 CALL iom_get( inum, 'time' , zdate_inc ) 341 342 CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb ) 342 343 CALL iom_get( inum, 'z_inc_datef', z_inc_datef ) … … 345 346 ! 346 347 IF(lwp) THEN 347 WRITE(numout,*) 348 WRITE(numout,*) 348 349 WRITE(numout,*) 'asm_inc_init : Assimilation increments valid between dates ', z_inc_dateb,' and ', z_inc_datef 349 350 WRITE(numout,*) '~~~~~~~~~~~~' … … 359 360 & ' not agree with Direct Initialization time' ) 360 361 361 IF ( ln_trainc ) THEN 362 IF ( ln_trainc ) THEN 362 363 CALL iom_get( inum, jpdom_auto, 'bckint', t_bkginc, 1 ) 363 364 CALL iom_get( inum, jpdom_auto, 'bckins', s_bkginc, 1 ) … … 371 372 ENDIF 372 373 373 IF ( ln_dyninc ) THEN 374 CALL iom_get( inum, jpdom_auto, 'bckinu', u_bkginc, 1 ) 375 CALL iom_get( inum, jpdom_auto, 'bckinv', v_bkginc, 1 ) 374 IF ( ln_dyninc ) THEN 375 CALL iom_get( inum, jpdom_auto, 'bckinu', u_bkginc, 1 ) 376 CALL iom_get( inum, jpdom_auto, 'bckinv', v_bkginc, 1 ) 376 377 ! Apply the masks 377 378 u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:) … … 382 383 WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0 383 384 ENDIF 384 385 385 386 IF ( ln_sshinc ) THEN 386 387 CALL iom_get( inum, jpdom_auto, 'bckineta', ssh_bkginc, 1 ) … … 408 409 IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN ! Apply divergence damping filter 409 410 ! !-------------------------------------- 410 ALLOCATE( zhdiv(jpi,jpj) ) 411 ALLOCATE( zhdiv(jpi,jpj) ) 411 412 ! 412 413 DO jt = 1, nn_divdmp … … 427 428 & + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 428 429 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) & 429 & + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 430 & + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 430 431 END_2D 431 432 END DO … … 433 434 END DO 434 435 ! 435 DEALLOCATE( zhdiv ) 436 DEALLOCATE( zhdiv ) 436 437 ! 437 438 ENDIF … … 454 455 CALL iom_open( c_asmdin, inum ) 455 456 ! 456 CALL iom_get( inum, 'rdastp', zdate_bkg ) 457 CALL iom_get( inum, 'rdastp', zdate_bkg ) 457 458 ! 458 459 IF(lwp) THEN 459 WRITE(numout,*) 460 WRITE(numout,*) 460 461 WRITE(numout,*) ' ==>>> Assimilation background state valid at : ', zdate_bkg 461 462 WRITE(numout,*) … … 466 467 & ' not agree with Direct Initialization time' ) 467 468 ! 468 IF ( ln_trainc ) THEN 469 IF ( ln_trainc ) THEN 469 470 CALL iom_get( inum, jpdom_auto, 'tn', t_bkg ) 470 471 CALL iom_get( inum, jpdom_auto, 'sn', s_bkg ) … … 473 474 ENDIF 474 475 ! 475 IF ( ln_dyninc ) THEN 476 IF ( ln_dyninc ) THEN 476 477 CALL iom_get( inum, jpdom_auto, 'un', u_bkg, cd_type = 'U', psgn = 1._wp ) 477 478 CALL iom_get( inum, jpdom_auto, 'vn', v_bkg, cd_type = 'V', psgn = 1._wp ) … … 501 502 ! 502 503 END SUBROUTINE asm_inc_init 503 504 504 505 505 506 SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs ) 506 507 !!---------------------------------------------------------------------- 507 508 !! *** ROUTINE tra_asm_inc *** 508 !! 509 !! 509 510 !! ** Purpose : Apply the tracer (T and S) assimilation increments 510 511 !! 511 512 !! ** Method : Direct initialization or Incremental Analysis Updating 512 513 !! 513 !! ** Action : 514 !! ** Action : 514 515 !!---------------------------------------------------------------------- 515 516 INTEGER , INTENT(in ) :: kt ! Current time step … … 518 519 ! 519 520 INTEGER :: ji, jj, jk 520 INTEGER :: it 521 INTEGER :: it, itile 521 522 REAL(wp) :: zincwgt ! IAU weight for current time step 522 REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values 523 !!---------------------------------------------------------------------- 524 ! 525 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 526 ! used to prevent the applied increments taking the temperature below the local freezing point 527 DO jk = 1, jpkm1 528 CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) ) 529 END DO 523 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: fzptnz ! 3d freezing point values 524 !!---------------------------------------------------------------------- 525 ! 526 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 527 ! used to prevent the applied increments taking the temperature below the local freezing point 528 IF( ln_temnofreeze ) THEN 529 DO jk = 1, jpkm1 530 CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) ) 531 END DO 532 ENDIF 530 533 ! 531 534 ! !-------------------------------------- … … 538 541 zincwgt = wgtiau(it) / rn_Dt ! IAU weight for the current time step 539 542 ! 540 IF(lwp) THEN 541 WRITE(numout,*) 542 WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 543 WRITE(numout,*) '~~~~~~~~~~~~' 543 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 544 IF(lwp) THEN 545 WRITE(numout,*) 546 WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 547 WRITE(numout,*) '~~~~~~~~~~~~' 548 ENDIF 544 549 ENDIF 545 550 ! … … 548 553 IF (ln_temnofreeze) THEN 549 554 ! Do not apply negative increments if the temperature will fall below freezing 550 WHERE(t_bkginc( :,:,jk) > 0.0_wp .OR. &551 & pts( :,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )552 pts( :,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt555 WHERE(t_bkginc(A2D(0),jk) > 0.0_wp .OR. & 556 & pts(A2D(0),jk,jp_tem,Kmm) + pts(A2D(0),jk,jp_tem,Krhs) + t_bkginc(A2D(0),jk) * wgtiau(it) > fzptnz(:,:,jk) ) 557 pts(A2D(0),jk,jp_tem,Krhs) = pts(A2D(0),jk,jp_tem,Krhs) + t_bkginc(A2D(0),jk) * zincwgt 553 558 END WHERE 554 559 ELSE 555 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 560 DO_2D( 0, 0, 0, 0 ) 561 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + t_bkginc(ji,jj,jk) * zincwgt 562 END_2D 556 563 ENDIF 557 564 IF (ln_salfix) THEN 558 565 ! Do not apply negative increments if the salinity will fall below a specified 559 566 ! minimum value salfixmin 560 WHERE(s_bkginc( :,:,jk) > 0.0_wp .OR. &561 & pts( :,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )562 pts( :,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt567 WHERE(s_bkginc(A2D(0),jk) > 0.0_wp .OR. & 568 & pts(A2D(0),jk,jp_sal,Kmm) + pts(A2D(0),jk,jp_sal,Krhs) + s_bkginc(A2D(0),jk) * wgtiau(it) > salfixmin ) 569 pts(A2D(0),jk,jp_sal,Krhs) = pts(A2D(0),jk,jp_sal,Krhs) + s_bkginc(A2D(0),jk) * zincwgt 563 570 END WHERE 564 571 ELSE 565 pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 572 DO_2D( 0, 0, 0, 0 ) 573 pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) + s_bkginc(ji,jj,jk) * zincwgt 574 END_2D 566 575 ENDIF 567 576 END DO … … 569 578 ENDIF 570 579 ! 571 IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work 572 DEALLOCATE( t_bkginc ) 573 DEALLOCATE( s_bkginc ) 580 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 581 IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work 582 DEALLOCATE( t_bkginc ) 583 DEALLOCATE( s_bkginc ) 584 ENDIF 574 585 ENDIF 575 586 ! !-------------------------------------- 576 587 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 577 588 ! !-------------------------------------- 578 ! 589 ! 579 590 IF ( kt == nitdin_r ) THEN 580 591 ! … … 584 595 IF (ln_temnofreeze) THEN 585 596 ! Do not apply negative increments if the temperature will fall below freezing 586 WHERE( t_bkginc( :,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) )587 pts( :,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)597 WHERE( t_bkginc(A2D(0),:) > 0.0_wp .OR. pts(A2D(0),:,jp_tem,Kmm) + t_bkginc(A2D(0),:) > fzptnz(:,:,:) ) 598 pts(A2D(0),:,jp_tem,Kmm) = t_bkg(A2D(0),:) + t_bkginc(A2D(0),:) 588 599 END WHERE 589 600 ELSE 590 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 601 DO_3D( 0, 0, 0, 0, 1, jpk ) 602 pts(ji,jj,jk,jp_tem,Kmm) = t_bkg(ji,jj,jk) + t_bkginc(ji,jj,jk) 603 END_3D 591 604 ENDIF 592 605 IF (ln_salfix) THEN 593 606 ! Do not apply negative increments if the salinity will fall below a specified 594 607 ! minimum value salfixmin 595 WHERE( s_bkginc( :,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin )596 pts( :,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)608 WHERE( s_bkginc(A2D(0),:) > 0.0_wp .OR. pts(A2D(0),:,jp_sal,Kmm) + s_bkginc(A2D(0),:) > salfixmin ) 609 pts(A2D(0),:,jp_sal,Kmm) = s_bkg(A2D(0),:) + s_bkginc(A2D(0),:) 597 610 END WHERE 598 611 ELSE 599 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 600 ENDIF 601 602 pts(:,:,:,:,Kbb) = pts(:,:,:,:,Kmm) ! Update before fields 612 DO_3D( 0, 0, 0, 0, 1, jpk ) 613 pts(ji,jj,jk,jp_sal,Kmm) = s_bkg(ji,jj,jk) + s_bkginc(ji,jj,jk) 614 END_3D 615 ENDIF 616 617 DO_3D( 0, 0, 0, 0, 1, jpk ) 618 pts(ji,jj,jk,:,Kbb) = pts(ji,jj,jk,:,Kmm) ! Update before fields 619 END_3D 603 620 604 621 CALL eos( pts(:,:,:,:,Kbb), rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities … … 607 624 !!gm 608 625 609 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 610 & CALL zps_hde ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient 611 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 612 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 613 & CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 614 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the last ocean level 615 616 DEALLOCATE( t_bkginc ) 617 DEALLOCATE( s_bkginc ) 618 DEALLOCATE( t_bkg ) 619 DEALLOCATE( s_bkg ) 620 ENDIF 621 ! 626 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from zps_hde*) 627 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 628 itile = ntile 629 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) ! Use full domain 630 631 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 632 & CALL zps_hde ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient 633 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 634 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 635 & CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 636 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the last ocean level 637 638 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = itile ) ! Revert to tile domain 639 ENDIF 640 641 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 642 DEALLOCATE( t_bkginc ) 643 DEALLOCATE( s_bkginc ) 644 DEALLOCATE( t_bkg ) 645 DEALLOCATE( s_bkg ) 646 ENDIF 647 ! 648 ENDIF 649 ! 622 650 ENDIF 623 651 ! Perhaps the following call should be in step … … 630 658 !!---------------------------------------------------------------------- 631 659 !! *** ROUTINE dyn_asm_inc *** 632 !! 660 !! 633 661 !! ** Purpose : Apply the dynamics (u and v) assimilation increments. 634 662 !! 635 663 !! ** Method : Direct initialization or Incremental Analysis Updating. 636 664 !! 637 !! ** Action : 665 !! ** Action : 638 666 !!---------------------------------------------------------------------- 639 667 INTEGER , INTENT( in ) :: kt ! ocean time-step index … … 656 684 ! 657 685 IF(lwp) THEN 658 WRITE(numout,*) 686 WRITE(numout,*) 659 687 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 660 688 WRITE(numout,*) '~~~~~~~~~~~~' … … 676 704 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 677 705 ! !----------------------------------------- 678 ! 706 ! 679 707 IF ( kt == nitdin_r ) THEN 680 708 ! … … 683 711 ! Initialize the now fields with the background + increment 684 712 puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:) 685 pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 713 pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 686 714 ! 687 715 puu(:,:,:,Kbb) = puu(:,:,:,Kmm) ! Update before fields … … 702 730 !!---------------------------------------------------------------------- 703 731 !! *** ROUTINE ssh_asm_inc *** 704 !! 732 !! 705 733 !! ** Purpose : Apply the sea surface height assimilation increment. 706 734 !! 707 735 !! ** Method : Direct initialization or Incremental Analysis Updating. 708 736 !! 709 !! ** Action : 737 !! ** Action : 710 738 !!---------------------------------------------------------------------- 711 739 INTEGER, INTENT(IN) :: kt ! Current time step … … 727 755 ! 728 756 IF(lwp) THEN 729 WRITE(numout,*) 757 WRITE(numout,*) 730 758 WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', & 731 759 & kt,' with IAU weight = ', wgtiau(it) … … 779 807 !! *** ROUTINE ssh_asm_div *** 780 808 !! 781 !! ** Purpose : ssh increment with z* is incorporated via a correction of the local divergence 809 !! ** Purpose : ssh increment with z* is incorporated via a correction of the local divergence 782 810 !! across all the water column 783 811 !! … … 795 823 REAL(wp), DIMENSION(:,:) , POINTER :: ztim ! local array 796 824 !!---------------------------------------------------------------------- 797 ! 825 ! 798 826 #if defined key_asminc 799 827 CALL ssh_asm_inc( kt, Kbb, Kmm ) !== (calculate increments) 800 828 ! 801 IF( ln_linssh ) THEN 829 IF( ln_linssh ) THEN 802 830 phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1) 803 ELSE 831 ELSE 804 832 ALLOCATE( ztim(jpi,jpj) ) 805 833 ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) ) 806 DO jk = 1, jpkm1 807 phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 834 DO jk = 1, jpkm1 835 phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 808 836 END DO 809 837 ! … … 818 846 !!---------------------------------------------------------------------- 819 847 !! *** ROUTINE seaice_asm_inc *** 820 !! 848 !! 821 849 !! ** Purpose : Apply the sea ice assimilation increment. 822 850 !! 823 851 !! ** Method : Direct initialization or Incremental Analysis Updating. 824 852 !! 825 !! ** Action : 853 !! ** Action : 826 854 !! 827 855 !!---------------------------------------------------------------------- … … 829 857 INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation 830 858 ! 859 INTEGER :: ji, jj 831 860 INTEGER :: it 832 861 REAL(wp) :: zincwgt ! IAU weight for current time step 833 862 #if defined key_si3 834 REAL(wp), DIMENSION( jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc863 REAL(wp), DIMENSION(A2D(nn_hls)) :: zofrld, zohicif, zseaicendg, zhicifinc 835 864 REAL(wp) :: zhicifmin = 0.5_wp ! ice minimum depth in metres 836 865 #endif … … 844 873 ! 845 874 it = kt - nit000 + 1 846 zincwgt = wgtiau(it) ! IAU weight for the current time step 875 zincwgt = wgtiau(it) ! IAU weight for the current time step 847 876 ! note this is not a tendency so should not be divided by rn_Dt (as with the tracer and other increments) 848 877 ! 849 IF(lwp) THEN 850 WRITE(numout,*) 851 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 852 WRITE(numout,*) '~~~~~~~~~~~~' 878 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 879 IF(lwp) THEN 880 WRITE(numout,*) 881 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 882 WRITE(numout,*) '~~~~~~~~~~~~' 883 ENDIF 853 884 ENDIF 854 885 ! … … 856 887 ! 857 888 #if defined key_si3 858 zofrld (:,:) = 1._wp - at_i(:,:) 859 zohicif(:,:) = hm_i(:,:) 860 ! 861 at_i (:,:) = 1. - MIN( MAX( 1.-at_i (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 862 at_i_b(:,:) = 1. - MIN( MAX( 1.-at_i_b(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 863 fr_i(:,:) = at_i(:,:) ! adjust ice fraction 864 ! 865 zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:)) ! find out actual sea ice nudge applied 889 DO_2D( 0, 0, 0, 0 ) 890 zofrld (ji,jj) = 1._wp - at_i(ji,jj) 891 zohicif(ji,jj) = hm_i(ji,jj) 892 ! 893 at_i (ji,jj) = 1. - MIN( MAX( 1.-at_i (ji,jj) - seaice_bkginc(ji,jj) * zincwgt, 0.0_wp), 1.0_wp) 894 at_i_b(ji,jj) = 1. - MIN( MAX( 1.-at_i_b(ji,jj) - seaice_bkginc(ji,jj) * zincwgt, 0.0_wp), 1.0_wp) 895 fr_i(ji,jj) = at_i(ji,jj) ! adjust ice fraction 896 ! 897 zseaicendg(ji,jj) = zofrld(ji,jj) - (1. - at_i(ji,jj)) ! find out actual sea ice nudge applied 898 END_2D 866 899 ! 867 900 ! Nudge sea ice depth to bring it up to a required minimum depth 868 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i( :,:) < zhicifmin )869 zhicifinc(:,:) = (zhicifmin - hm_i( :,:)) * zincwgt901 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(A2D(0)) < zhicifmin ) 902 zhicifinc(:,:) = (zhicifmin - hm_i(A2D(0))) * zincwgt 870 903 ELSEWHERE 871 904 zhicifinc(:,:) = 0.0_wp … … 873 906 ! 874 907 ! nudge ice depth 875 hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:) 908 DO_2D( 0, 0, 0, 0 ) 909 hm_i (ji,jj) = hm_i (ji,jj) + zhicifinc(ji,jj) 910 END_2D 876 911 ! 877 912 ! seaice salinity balancing (to add) … … 880 915 #if defined key_cice && defined key_asminc 881 916 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 882 ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rn_Dt 883 #endif 884 ! 885 IF ( kt == nitiaufin_r ) THEN 886 DEALLOCATE( seaice_bkginc ) 917 DO_2D( 0, 0, 0, 0 ) 918 ndaice_da(ji,jj) = seaice_bkginc(ji,jj) * zincwgt / rn_Dt 919 END_2D 920 #endif 921 ! 922 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 923 IF ( kt == nitiaufin_r ) THEN 924 DEALLOCATE( seaice_bkginc ) 925 ENDIF 887 926 ENDIF 888 927 ! … … 890 929 ! 891 930 #if defined key_cice && defined key_asminc 892 ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 931 DO_2D( 0, 0, 0, 0 ) 932 ndaice_da(ji,jj) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 933 END_2D 893 934 #endif 894 935 ! … … 905 946 ! 906 947 #if defined key_si3 907 zofrld (:,:) = 1._wp - at_i(:,:) 908 zohicif(:,:) = hm_i(:,:) 909 ! 910 ! Initialize the now fields the background + increment 911 at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 912 at_i_b(:,:) = at_i(:,:) 913 fr_i(:,:) = at_i(:,:) ! adjust ice fraction 914 ! 915 zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:)) ! find out actual sea ice nudge applied 948 DO_2D( 0, 0, 0, 0 ) 949 zofrld (ji,jj) = 1._wp - at_i(ji,jj) 950 zohicif(ji,jj) = hm_i(ji,jj) 951 ! 952 ! Initialize the now fields the background + increment 953 at_i(ji,jj) = 1. - MIN( MAX( 1.-at_i(ji,jj) - seaice_bkginc(ji,jj), 0.0_wp), 1.0_wp) 954 at_i_b(ji,jj) = at_i(ji,jj) 955 fr_i(ji,jj) = at_i(ji,jj) ! adjust ice fraction 956 ! 957 zseaicendg(ji,jj) = zofrld(ji,jj) - (1. - at_i(ji,jj)) ! find out actual sea ice nudge applied 958 END_2D 916 959 ! 917 960 ! Nudge sea ice depth to bring it up to a required minimum depth 918 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i( :,:) < zhicifmin )919 zhicifinc(:,:) = zhicifmin - hm_i( :,:)961 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(A2D(0)) < zhicifmin ) 962 zhicifinc(:,:) = zhicifmin - hm_i(A2D(0)) 920 963 ELSEWHERE 921 964 zhicifinc(:,:) = 0.0_wp … … 923 966 ! 924 967 ! nudge ice depth 925 hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:) 968 DO_2D( 0, 0, 0, 0 ) 969 hm_i(ji,jj) = hm_i (ji,jj) + zhicifinc(ji,jj) 970 END_2D 926 971 ! 927 972 ! seaice salinity balancing (to add) … … 930 975 #if defined key_cice && defined key_asminc 931 976 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 932 ndaice_da(:,:) = seaice_bkginc(:,:) / rn_Dt 933 #endif 934 IF ( .NOT. PRESENT(kindic) ) THEN 935 DEALLOCATE( seaice_bkginc ) 936 END IF 977 DO_2D( 0, 0, 0, 0 ) 978 ndaice_da(ji,jj) = seaice_bkginc(ji,jj) / rn_Dt 979 END_2D 980 #endif 981 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 982 IF ( .NOT. PRESENT(kindic) ) THEN 983 DEALLOCATE( seaice_bkginc ) 984 END IF 985 ENDIF 937 986 ! 938 987 ELSE 939 988 ! 940 989 #if defined key_cice && defined key_asminc 941 ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 990 DO_2D( 0, 0, 0, 0 ) 991 ndaice_da(ji,jj) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 992 END_2D 942 993 #endif 943 994 ! … … 946 997 !#if defined defined key_si3 || defined key_cice 947 998 ! 948 ! IF (ln_seaicebal ) THEN 999 ! IF (ln_seaicebal ) THEN 949 1000 ! !! balancing salinity increments 950 1001 ! !! simple case from limflx.F90 (doesn't include a mass flux) … … 958 1009 ! 959 1010 ! DO jj = 1, jpj 960 ! DO ji = 1, jpi 1011 ! DO ji = 1, jpi 961 1012 ! ! calculate change in ice and snow mass per unit area 962 1013 ! ! positive values imply adding salt to the ocean (results from ice formation) … … 969 1020 ! 970 1021 ! ! prevent small mld 971 ! ! less than 10m can cause salinity instability 1022 ! ! less than 10m can cause salinity instability 972 1023 ! IF (mld < 10) mld=10 973 1024 ! 974 ! ! set to bottom of a level 1025 ! ! set to bottom of a level 975 1026 ! DO jk = jpk-1, 2, -1 976 1027 ! IF ((mld > gdepw(ji,jj,jk,Kmm)) .and. (mld < gdepw(ji,jj,jk+1,Kmm))) THEN … … 981 1032 ! 982 1033 ! ! avoid applying salinity balancing in shallow water or on land 983 ! ! 1034 ! ! 984 1035 ! 985 1036 ! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) … … 992 1043 ! 993 1044 ! ! put increments in for levels in the mixed layer 994 ! ! but prevent salinity below a threshold value 995 ! 996 ! DO jk = 1, jkmax 997 ! 998 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 1045 ! ! but prevent salinity below a threshold value 1046 ! 1047 ! DO jk = 1, jkmax 1048 ! 1049 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 999 1050 ! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 1000 1051 ! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn … … 1007 1058 ! ! 1008 1059 ! !! Adjust fsalt. A +ve fsalt means adding salt to ocean 1009 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1010 ! !! 1011 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1060 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1061 ! !! 1062 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1012 1063 ! !! ! E-P (kg m-2 s-2) 1013 1064 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2) … … 1022 1073 ! 1023 1074 END SUBROUTINE seaice_asm_inc 1024 1075 1025 1076 !!====================================================================== 1026 1077 END MODULE asminc
Note: See TracChangeset
for help on using the changeset viewer.
