- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/ASM/asminc.F90 (modified) (31 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ASM/asminc.F90
r12178 r12928 94 94 95 95 !! * Substitutions 96 # include " vectopt_loop_substitute.h90"96 # include "do_loop_substitute.h90" 97 97 !!---------------------------------------------------------------------- 98 98 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 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 … … 145 147 ln_temnofreeze = .FALSE. 146 148 147 REWIND( numnam_ref ) ! Namelist nam_asminc in reference namelist : Assimilation increment148 149 READ ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901) 149 150 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in reference namelist' ) 150 REWIND( numnam_cfg ) ! Namelist nam_asminc in configuration namelist : Assimilation increment151 151 READ ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 ) 152 152 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_asminc in configuration namelist' ) … … 413 413 DO jk = 1, jpkm1 ! zhdiv = e1e1 * div 414 414 zhdiv(:,:) = 0._wp 415 DO jj = 2, jpjm1 416 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) 421 END DO 422 END DO 415 DO_2D_00_00 416 zhdiv(ji,jj) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * u_bkginc(ji ,jj,jk) & 417 & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk) & 418 & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * v_bkginc(ji,jj ,jk) & 419 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk) ) / e3t(ji,jj,jk,Kmm) 420 END_2D 423 421 CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. ) ! lateral boundary cond. (no sign change) 424 422 ! 425 DO jj = 2, jpjm1 426 DO ji = fs_2, fs_jpim1 ! vector opt. 427 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) & 428 & + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 429 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) & 430 & + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 431 END DO 432 END DO 423 DO_2D_00_00 424 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) & 425 & + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 426 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) & 427 & + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 428 END_2D 433 429 END DO 434 430 ! … … 491 487 ENDIF 492 488 ! 493 IF(lwp) WRITE(numout,*) ' ==>>> Euler time step switch is ', neuler489 IF(lwp) WRITE(numout,*) ' ==>>> Euler time step switch is ', l_1st_euler 494 490 ! 495 491 IF( lk_asminc ) THEN !== data assimilation ==! 496 IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 ) ! Output background fields492 IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1, Kmm ) ! Output background fields 497 493 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 )! SSH494 IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1, Kbb, Kmm, ts , Krhs ) ! Tracers 495 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1, Kbb, Kmm, uu, vv, Krhs ) ! Dynamics 496 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1, Kbb, Kmm ) ! SSH 501 497 ENDIF 502 498 ENDIF … … 505 501 506 502 507 SUBROUTINE tra_asm_inc( kt )503 SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs ) 508 504 !!---------------------------------------------------------------------- 509 505 !! *** ROUTINE tra_asm_inc *** … … 515 511 !! ** Action : 516 512 !!---------------------------------------------------------------------- 517 INTEGER, INTENT(IN) :: kt ! Current time step 513 INTEGER , INTENT(in ) :: kt ! Current time step 514 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! Time level indices 515 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 518 516 ! 519 517 INTEGER :: ji, jj, jk … … 526 524 ! used to prevent the applied increments taking the temperature below the local freezing point 527 525 DO jk = 1, jpkm1 528 CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), gdept_n(:,:,jk) )526 CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) ) 529 527 END DO 530 528 ! … … 536 534 ! 537 535 it = kt - nit000 + 1 538 zincwgt = wgtiau(it) / r dt ! IAU weight for the current time step536 zincwgt = wgtiau(it) / rn_Dt ! IAU weight for the current time step 539 537 ! 540 538 IF(lwp) THEN … … 549 547 ! Do not apply negative increments if the temperature will fall below freezing 550 548 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) * zincwgt549 & pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 550 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 553 551 END WHERE 554 552 ELSE 555 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt553 pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 556 554 ENDIF 557 555 IF (ln_salfix) THEN … … 559 557 ! minimum value salfixmin 560 558 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) * zincwgt559 & pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 560 pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 563 561 END WHERE 564 562 ELSE 565 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt563 pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 566 564 ENDIF 567 565 END DO … … 579 577 IF ( kt == nitdin_r ) THEN 580 578 ! 581 neuler = 0! Force Euler forward step579 l_1st_euler = .TRUE. ! Force Euler forward step 582 580 ! 583 581 ! Initialize the now fields with the background + increment 584 582 IF (ln_temnofreeze) THEN 585 583 ! 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(:,:,:)584 WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 585 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 588 586 END WHERE 589 587 ELSE 590 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)588 pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 591 589 ENDIF 592 590 IF (ln_salfix) THEN 593 591 ! Do not apply negative increments if the salinity will fall below a specified 594 592 ! 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(:,:,:)593 WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin ) 594 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 597 595 END WHERE 598 596 ELSE 599 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)597 pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 600 598 ENDIF 601 599 602 tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields603 604 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities600 pts(:,:,:,:,Kbb) = pts(:,:,:,:,Kmm) ! Update before fields 601 602 CALL eos( pts(:,:,:,:,Kbb), rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities 605 603 !!gm fabien 606 ! CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities604 ! CALL eos( pts(:,:,:,:,Kbb), rhd, rhop ) ! Before potential and in situ densities 607 605 !!gm 608 606 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 level607 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 608 & CALL zps_hde ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient 609 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 610 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 611 & CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 612 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the last ocean level 615 613 616 614 DEALLOCATE( t_bkginc ) … … 627 625 628 626 629 SUBROUTINE dyn_asm_inc( kt )627 SUBROUTINE dyn_asm_inc( kt, Kbb, Kmm, puu, pvv, Krhs ) 630 628 !!---------------------------------------------------------------------- 631 629 !! *** ROUTINE dyn_asm_inc *** … … 637 635 !! ** Action : 638 636 !!---------------------------------------------------------------------- 639 INTEGER, INTENT(IN) :: kt ! Current time step 637 INTEGER , INTENT( in ) :: kt ! ocean time-step index 638 INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices 639 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 640 640 ! 641 641 INTEGER :: jk … … 651 651 ! 652 652 it = kt - nit000 + 1 653 zincwgt = wgtiau(it) / r dt ! IAU weight for the current time step653 zincwgt = wgtiau(it) / rn_Dt ! IAU weight for the current time step 654 654 ! 655 655 IF(lwp) THEN … … 661 661 ! Update the dynamic tendencies 662 662 DO jk = 1, jpkm1 663 ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt664 va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt663 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + u_bkginc(:,:,jk) * zincwgt 664 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + v_bkginc(:,:,jk) * zincwgt 665 665 END DO 666 666 ! … … 677 677 IF ( kt == nitdin_r ) THEN 678 678 ! 679 neuler = 0! Force Euler forward step679 l_1st_euler = .TRUE. ! Force Euler forward step 680 680 ! 681 681 ! 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(:,:,:)682 puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:) 683 pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 684 ! 685 puu(:,:,:,Kbb) = puu(:,:,:,Kmm) ! Update before fields 686 pvv(:,:,:,Kbb) = pvv(:,:,:,Kmm) 687 687 ! 688 688 DEALLOCATE( u_bkg ) … … 697 697 698 698 699 SUBROUTINE ssh_asm_inc( kt )699 SUBROUTINE ssh_asm_inc( kt, Kbb, Kmm ) 700 700 !!---------------------------------------------------------------------- 701 701 !! *** ROUTINE ssh_asm_inc *** … … 707 707 !! ** Action : 708 708 !!---------------------------------------------------------------------- 709 INTEGER, INTENT(IN) :: kt ! Current time step 709 INTEGER, INTENT(IN) :: kt ! Current time step 710 INTEGER, INTENT(IN) :: Kbb, Kmm ! Current time step 710 711 ! 711 712 INTEGER :: it … … 721 722 ! 722 723 it = kt - nit000 + 1 723 zincwgt = wgtiau(it) / r dt ! IAU weight for the current time step724 zincwgt = wgtiau(it) / rn_Dt ! IAU weight for the current time step 724 725 ! 725 726 IF(lwp) THEN … … 752 753 IF ( kt == nitdin_r ) THEN 753 754 ! 754 neuler = 0! Force Euler forward step755 ! 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????755 l_1st_euler = .TRUE. ! Force Euler forward step 756 ! 757 ssh(:,:,Kmm) = ssh_bkg(:,:) + ssh_bkginc(:,:) ! Initialize the now fields the background + increment 758 ! 759 ssh(:,:,Kbb) = ssh(:,:,Kmm) ! Update before fields 760 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 761 !!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ???? 761 762 ! 762 763 DEALLOCATE( ssh_bkg ) … … 770 771 771 772 772 SUBROUTINE ssh_asm_div( kt, phdivn )773 SUBROUTINE ssh_asm_div( kt, Kbb, Kmm, phdivn ) 773 774 !!---------------------------------------------------------------------- 774 775 !! *** ROUTINE ssh_asm_div *** … … 784 785 !!---------------------------------------------------------------------- 785 786 INTEGER, INTENT(IN) :: kt ! ocean time-step index 787 INTEGER, INTENT(IN) :: Kbb, Kmm ! time level indices 786 788 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence 787 789 !! … … 791 793 ! 792 794 #if defined key_asminc 793 CALL ssh_asm_inc( kt ) !== (calculate increments)795 CALL ssh_asm_inc( kt, Kbb, Kmm ) !== (calculate increments) 794 796 ! 795 797 IF( ln_linssh ) THEN 796 phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t _n(:,:,1) * tmask(:,:,1)798 phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1) 797 799 ELSE 798 800 ALLOCATE( ztim(jpi,jpj) ) 799 ztim(:,:) = ssh_iau(:,:) / ( ht _n(:,:) + 1.0 - ssmask(:,:) )801 ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) ) 800 802 DO jk = 1, jpkm1 801 803 phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) … … 839 841 it = kt - nit000 + 1 840 842 zincwgt = wgtiau(it) ! IAU weight for the current time step 841 ! note this is not a tendency so should not be divided by r dt (as with the tracer and other increments)843 ! note this is not a tendency so should not be divided by rn_Dt (as with the tracer and other increments) 842 844 ! 843 845 IF(lwp) THEN … … 874 876 #if defined key_cice && defined key_asminc 875 877 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 876 ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / r dt878 ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rn_Dt 877 879 #endif 878 880 ! … … 894 896 IF ( kt == nitdin_r ) THEN 895 897 ! 896 neuler = 0! Force Euler forward step898 l_1st_euler = .TRUE. ! Force Euler forward step 897 899 ! 898 900 ! Sea-ice : SI3 case … … 924 926 #if defined key_cice && defined key_asminc 925 927 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 926 ndaice_da(:,:) = seaice_bkginc(:,:) / r dt928 ndaice_da(:,:) = seaice_bkginc(:,:) / rn_Dt 927 929 #endif 928 930 IF ( .NOT. PRESENT(kindic) ) THEN … … 957 959 ! ! fwf : ice formation and melting 958 960 ! 959 ! zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*r dt961 ! zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rn_Dt 960 962 ! 961 963 ! ! change salinity down to mixed layer depth … … 998 1000 ! 999 1001 ! ! ! salt exchanges at the ice/ocean interface 1000 ! ! zpmess = zfons / r dt_ice ! rdt_ice is ice timestep1002 ! ! zpmess = zfons / rDt_ice ! rDt_ice is ice timestep 1001 1003 ! ! 1002 1004 ! !! Adjust fsalt. A +ve fsalt means adding salt to ocean
Note: See TracChangeset
for help on using the changeset viewer.
