- Timestamp:
- 2019-11-22T15:29:17+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
OCE/ASM/asminc.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Property svn:mergeinfo deleted
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ASM/asminc.F90
r11536 r11949 102 102 CONTAINS 103 103 104 SUBROUTINE asm_inc_init 104 SUBROUTINE asm_inc_init( Kbb, Kmm, Krhs ) 105 105 !!---------------------------------------------------------------------- 106 106 !! *** ROUTINE asm_inc_init *** … … 112 112 !! ** Action : 113 113 !!---------------------------------------------------------------------- 114 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices 115 ! 114 116 INTEGER :: ji, jj, jk, jt ! dummy loop indices 115 117 INTEGER :: imid, inum ! local integers … … 415 417 DO jj = 2, jpjm1 416 418 DO ji = fs_2, fs_jpim1 ! vector opt. 417 zhdiv(ji,jj) = ( e2u(ji ,jj) * e3u _n(ji ,jj,jk) * u_bkginc(ji ,jj,jk) &418 & - e2u(ji-1,jj) * e3u _n(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk) &419 & + e1v(ji,jj ) * e3v _n(ji,jj ,jk) * v_bkginc(ji,jj ,jk) &420 & - e1v(ji,jj-1) * e3v _n(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk) ) / e3t_n(ji,jj,jk)419 zhdiv(ji,jj) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * u_bkginc(ji ,jj,jk) & 420 & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk) & 421 & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * v_bkginc(ji,jj ,jk) & 422 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk) ) / e3t(ji,jj,jk,Kmm) 421 423 END DO 422 424 END DO … … 494 496 ! 495 497 IF( lk_asminc ) THEN !== data assimilation ==! 496 IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 ) ! Output background fields498 IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1, Kmm ) ! Output background fields 497 499 IF( ln_asmdin ) THEN ! Direct initialization 498 IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 ) ! Tracers499 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 ) ! Dynamics500 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )! SSH500 IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1, Kbb, Kmm, ts , Krhs ) ! Tracers 501 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1, Kbb, Kmm, uu, vv, Krhs ) ! Dynamics 502 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1, Kbb, Kmm ) ! SSH 501 503 ENDIF 502 504 ENDIF … … 505 507 506 508 507 SUBROUTINE tra_asm_inc( kt )509 SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs ) 508 510 !!---------------------------------------------------------------------- 509 511 !! *** ROUTINE tra_asm_inc *** … … 515 517 !! ** Action : 516 518 !!---------------------------------------------------------------------- 517 INTEGER, INTENT(IN) :: kt ! Current time step 519 INTEGER , INTENT(in ) :: kt ! Current time step 520 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! Time level indices 521 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 518 522 ! 519 523 INTEGER :: ji, jj, jk … … 526 530 ! used to prevent the applied increments taking the temperature below the local freezing point 527 531 DO jk = 1, jpkm1 528 CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), gdept_n(:,:,jk) )532 CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) ) 529 533 END DO 530 534 ! … … 549 553 ! Do not apply negative increments if the temperature will fall below freezing 550 554 WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. & 551 & tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )552 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt555 & pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 556 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 553 557 END WHERE 554 558 ELSE 555 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt559 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 556 560 ENDIF 557 561 IF (ln_salfix) THEN … … 559 563 ! minimum value salfixmin 560 564 WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. & 561 & tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )562 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt565 & pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 566 pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 563 567 END WHERE 564 568 ELSE 565 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt569 pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 566 570 ENDIF 567 571 END DO … … 584 588 IF (ln_temnofreeze) THEN 585 589 ! Do not apply negative increments if the temperature will fall below freezing 586 WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )587 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)590 WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 591 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 588 592 END WHERE 589 593 ELSE 590 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)594 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 591 595 ENDIF 592 596 IF (ln_salfix) THEN 593 597 ! Do not apply negative increments if the salinity will fall below a specified 594 598 ! minimum value salfixmin 595 WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )596 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)599 WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin ) 600 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 597 601 END WHERE 598 602 ELSE 599 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)603 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 600 604 ENDIF 601 605 602 tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields603 604 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities606 pts(:,:,:,:,Kbb) = pts(:,:,:,:,Kmm) ! Update before fields 607 608 CALL eos( pts(:,:,:,:,Kbb), rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities 605 609 !!gm fabien 606 ! CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities610 ! CALL eos( pts(:,:,:,:,Kbb), rhd, rhop ) ! Before potential and in situ densities 607 611 !!gm 608 612 609 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) &610 & CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient611 & rhd, gru , grv ) ! of t, s, rd at the last ocean level612 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) &613 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi, &! Partial steps for top cell (ISF)614 & rhd, gru , grv , grui, grvi )! of t, s, rd at the last ocean level613 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 614 & CALL zps_hde ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient 615 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 616 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 617 & CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 618 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the last ocean level 615 619 616 620 DEALLOCATE( t_bkginc ) … … 627 631 628 632 629 SUBROUTINE dyn_asm_inc( kt )633 SUBROUTINE dyn_asm_inc( kt, Kbb, Kmm, puu, pvv, Krhs ) 630 634 !!---------------------------------------------------------------------- 631 635 !! *** ROUTINE dyn_asm_inc *** … … 637 641 !! ** Action : 638 642 !!---------------------------------------------------------------------- 639 INTEGER, INTENT(IN) :: kt ! Current time step 643 INTEGER , INTENT( in ) :: kt ! ocean time-step index 644 INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices 645 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 640 646 ! 641 647 INTEGER :: jk … … 661 667 ! Update the dynamic tendencies 662 668 DO jk = 1, jpkm1 663 ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt664 va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt669 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + u_bkginc(:,:,jk) * zincwgt 670 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + v_bkginc(:,:,jk) * zincwgt 665 671 END DO 666 672 ! … … 680 686 ! 681 687 ! Initialize the now fields with the background + increment 682 un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)683 vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)684 ! 685 ub(:,:,:) = un(:,:,:) ! Update before fields686 vb(:,:,:) = vn(:,:,:)688 puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:) 689 pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 690 ! 691 puu(:,:,:,Kbb) = puu(:,:,:,Kmm) ! Update before fields 692 pvv(:,:,:,Kbb) = pvv(:,:,:,Kmm) 687 693 ! 688 694 DEALLOCATE( u_bkg ) … … 697 703 698 704 699 SUBROUTINE ssh_asm_inc( kt )705 SUBROUTINE ssh_asm_inc( kt, Kbb, Kmm ) 700 706 !!---------------------------------------------------------------------- 701 707 !! *** ROUTINE ssh_asm_inc *** … … 707 713 !! ** Action : 708 714 !!---------------------------------------------------------------------- 709 INTEGER, INTENT(IN) :: kt ! Current time step 715 INTEGER, INTENT(IN) :: kt ! Current time step 716 INTEGER, INTENT(IN) :: Kbb, Kmm ! Current time step 710 717 ! 711 718 INTEGER :: it … … 754 761 neuler = 0 ! Force Euler forward step 755 762 ! 756 ssh n(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) ! Initialize the now fields the background + increment757 ! 758 ssh b(:,:) = sshn(:,:) ! Update before fields759 e3t _b(:,:,:) = e3t_n(:,:,:)760 !!gm why not e3u _b, e3v_b, gdept_b????763 ssh(:,:,Kmm) = ssh_bkg(:,:) + ssh_bkginc(:,:) ! Initialize the now fields the background + increment 764 ! 765 ssh(:,:,Kbb) = ssh(:,:,Kmm) ! Update before fields 766 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 767 !!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ???? 761 768 ! 762 769 DEALLOCATE( ssh_bkg ) … … 770 777 771 778 772 SUBROUTINE ssh_asm_div( kt, phdivn )779 SUBROUTINE ssh_asm_div( kt, Kbb, Kmm, phdivn ) 773 780 !!---------------------------------------------------------------------- 774 781 !! *** ROUTINE ssh_asm_div *** … … 784 791 !!---------------------------------------------------------------------- 785 792 INTEGER, INTENT(IN) :: kt ! ocean time-step index 793 INTEGER, INTENT(IN) :: Kbb, Kmm ! time level indices 786 794 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence 787 795 !! … … 791 799 ! 792 800 #if defined key_asminc 793 CALL ssh_asm_inc( kt ) !== (calculate increments)801 CALL ssh_asm_inc( kt, Kbb, Kmm ) !== (calculate increments) 794 802 ! 795 803 IF( ln_linssh ) THEN 796 phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t _n(:,:,1) * tmask(:,:,1)804 phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1) 797 805 ELSE 798 806 ALLOCATE( ztim(jpi,jpj) ) 799 ztim(:,:) = ssh_iau(:,:) / ( ht _n(:,:) + 1.0 - ssmask(:,:) )807 ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) ) 800 808 DO jk = 1, jpkm1 801 809 phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk)
Note: See TracChangeset
for help on using the changeset viewer.
