Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/DIA
- Timestamp:
- 2015-12-21T12:35:23+01:00 (8 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/DIA
- Files:
-
- 2 deleted
- 9 edited
- 3 copied
-
dia25h.F90 (copied) (copied from branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90)
-
diaar5.F90 (modified) (7 diffs)
-
diacfl.F90 (copied) (copied from branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90)
-
diadct.F90 (modified) (19 diffs)
-
diadimg.F90 (deleted)
-
diafwb.F90 (modified) (11 diffs)
-
diaharm.F90 (modified) (8 diffs)
-
diahsb.F90 (modified) (13 diffs)
-
diahth.F90 (modified) (11 diffs)
-
dianam.F90 (modified) (2 diffs)
-
diaptr.F90 (modified) (2 diffs)
-
diatmb.F90 (copied) (copied from branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90)
-
diawri.F90 (modified) (35 diffs)
-
diawri_dimg.h90 (deleted)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
r5836 r6140 40 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sn0 ! initial salinity 41 41 42 !! * Substitutions43 # include "domzgr_substitute.h90"44 42 !!---------------------------------------------------------------------- 45 43 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 99 97 ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) ! thermosteric ssh 100 98 ztsn(:,:,:,jp_sal) = sn0(:,:,:) 101 CALL eos( ztsn, zrhd, fsdept_n(:,:,:) ) ! now in situ density using initial salinity99 CALL eos( ztsn, zrhd, gdept_n(:,:,:) ) ! now in situ density using initial salinity 102 100 ! 103 101 zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice 104 102 DO jk = 1, jpkm1 105 zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)106 END DO 107 IF( .NOT.lk_vvl) THEN108 IF 103 zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk) 104 END DO 105 IF( ln_linssh ) THEN 106 IF( ln_isfcav ) THEN 109 107 DO ji=1,jpi 110 108 DO jj=1,jpj … … 115 113 zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 116 114 END IF 115 !!gm 116 !!gm riceload should be added in both ln_linssh=T or F, no? 117 !!gm 117 118 END IF 118 119 ! … … 123 124 124 125 ! ! steric sea surface height 125 CALL eos( tsn, zrhd, zrhop, fsdept_n(:,:,:) ) ! now in situ and potential density126 CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) ) ! now in situ and potential density 126 127 zrhop(:,:,jpk) = 0._wp 127 128 CALL iom_put( 'rhop', zrhop ) … … 129 130 zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice 130 131 DO jk = 1, jpkm1 131 zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)132 END DO 133 IF( .NOT.lk_vvl) THEN132 zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk) 133 END DO 134 IF( ln_linssh ) THEN 134 135 IF ( ln_isfcav ) THEN 135 136 DO ji=1,jpi … … 159 160 DO jj = 1, jpj 160 161 DO ji = 1, jpi 161 zztmp = area(ji,jj) * fse3t(ji,jj,jk)162 zztmp = area(ji,jj) * e3t_n(ji,jj,jk) 162 163 ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem) 163 164 zsal = zsal + zztmp * tsn(ji,jj,jk,jp_sal) … … 165 166 END DO 166 167 END DO 167 IF( .NOT.lk_vvl) THEN168 IF 168 IF( ln_linssh ) THEN 169 IF( ln_isfcav ) THEN 169 170 DO ji=1,jpi 170 171 DO jj=1,jpj -
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
r5505 r6140 1 1 MODULE diadct 2 !!===================================================================== 3 !! *** MODULE diadct *** 4 !! Ocean diagnostics: Compute the transport trough a sec. 5 !!=============================================================== 6 !! History : 7 !! 8 !! original : 02/99 (Y Drillet) 9 !! addition : 10/01 (Y Drillet, R Bourdalle Badie) 10 !! : 10/05 (M Laborie) F90 11 !! addition : 04/07 (G Garric) Ice sections 12 !! bugfix : 04/07 (C Bricaud) test on sec%nb_point 13 !! initialisation of ztransp1,ztransp2,... 14 !! nemo_v_3_4: 09/2011 (C Bricaud) 15 !! 16 !! 17 !!---------------------------------------------------------------------- 2 !!====================================================================== 3 !! *** MODULE diadct *** 4 !! Ocean diagnostics: Compute the transport trough a sec. 5 !!====================================================================== 6 !! History : OPA ! 02/1999 (Y Drillet) original code 7 !! ! 10/2001 (Y Drillet, R Bourdalle Badie) 8 !! NEMO 1.0 ! 10/2005 (M Laborie) F90 9 !! 3.0 ! 04/2007 (G Garric) Ice sections 10 !! - ! 04/2007 (C Bricaud) test on sec%nb_point, initialisation of ztransp1,ztransp2,... 11 !! 3.4 ! 09/2011 (C Bricaud) 12 !!---------------------------------------------------------------------- 18 13 #if defined key_diadct 19 !!---------------------------------------------------------------------- 20 !! 'key_diadct' : 21 !!---------------------------------------------------------------------- 22 !!---------------------------------------------------------------------- 23 !! dia_dct : Compute the transport through a sec. 24 !! dia_dct_init : Read namelist. 25 !! readsec : Read sections description and pathway 26 !! removepoints : Remove points which are common to 2 procs 27 !! transport : Compute transport for each sections 28 !! dia_dct_wri : Write tranports results in ascii files 29 !! interp : Compute temperature/salinity/density at U-point or V-point 30 !! 31 !!---------------------------------------------------------------------- 32 !! * Modules used 33 USE oce ! ocean dynamics and tracers 34 USE dom_oce ! ocean space and time domain 35 USE phycst ! physical constants 36 USE in_out_manager ! I/O manager 37 USE daymod ! calendar 38 USE dianam ! build name of file 39 USE lib_mpp ! distributed memory computing library 14 !!---------------------------------------------------------------------- 15 !! 'key_diadct' : 16 !!---------------------------------------------------------------------- 17 !!---------------------------------------------------------------------- 18 !! dia_dct : Compute the transport through a sec. 19 !! dia_dct_init : Read namelist. 20 !! readsec : Read sections description and pathway 21 !! removepoints : Remove points which are common to 2 procs 22 !! transport : Compute transport for each sections 23 !! dia_dct_wri : Write tranports results in ascii files 24 !! interp : Compute temperature/salinity/density at U-point or V-point 25 !! 26 !!---------------------------------------------------------------------- 27 USE oce ! ocean dynamics and tracers 28 USE dom_oce ! ocean space and time domain 29 USE phycst ! physical constants 30 USE in_out_manager ! I/O manager 31 USE daymod ! calendar 32 USE dianam ! build name of file 33 USE lib_mpp ! distributed memory computing library 40 34 #if defined key_lim2 41 USE ice_235 USE ice_2 42 36 #endif 43 37 #if defined key_lim3 44 USE ice38 USE ice 45 39 #endif 46 USE domvvl 47 USE timing ! preformance summary 48 USE wrk_nemo ! working arrays 49 50 IMPLICIT NONE 51 PRIVATE 52 53 !! * Routine accessibility 54 PUBLIC dia_dct ! routine called by step.F90 55 PUBLIC dia_dct_init ! routine called by opa.F90 56 PUBLIC diadct_alloc ! routine called by nemo_init in nemogcm.F90 57 PRIVATE readsec 58 PRIVATE removepoints 59 PRIVATE transport 60 PRIVATE dia_dct_wri 61 62 #include "domzgr_substitute.h90" 63 64 !! * Shared module variables 65 LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .TRUE. !: model-data diagnostics flag 66 67 !! * Module variables 68 INTEGER :: nn_dct ! Frequency of computation 69 INTEGER :: nn_dctwri ! Frequency of output 70 INTEGER :: nn_secdebug ! Number of the section to debug 40 USE domvvl 41 USE timing ! preformance summary 42 USE wrk_nemo ! working arrays 43 44 IMPLICIT NONE 45 PRIVATE 46 47 PUBLIC dia_dct ! routine called by step.F90 48 PUBLIC dia_dct_init ! routine called by opa.F90 49 PUBLIC diadct_alloc ! routine called by nemo_init in nemogcm.F90 50 PRIVATE readsec 51 PRIVATE removepoints 52 PRIVATE transport 53 PRIVATE dia_dct_wri 54 55 LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .TRUE. !: model-data diagnostics flag 56 57 INTEGER :: nn_dct ! Frequency of computation 58 INTEGER :: nn_dctwri ! Frequency of output 59 INTEGER :: nn_secdebug ! Number of the section to debug 71 60 72 INTEGER, PARAMETER :: nb_class_max = 1073 INTEGER, PARAMETER :: nb_sec_max = 15074 INTEGER, PARAMETER :: nb_point_max = 200075 INTEGER, PARAMETER :: nb_type_class = 1076 INTEGER, PARAMETER :: nb_3d_vars = 377 INTEGER, PARAMETER :: nb_2d_vars = 278 INTEGER :: nb_sec79 80 TYPE POINT_SECTION81 INTEGER :: I,J82 END TYPE POINT_SECTION83 84 TYPE COORD_SECTION85 REAL(wp) :: lon,lat86 END TYPE COORD_SECTION87 88 TYPE SECTION89 CHARACTER(len=60) :: name ! name of the sec90 LOGICAL :: llstrpond ! true if you want the computation of salt and61 INTEGER, PARAMETER :: nb_class_max = 10 62 INTEGER, PARAMETER :: nb_sec_max = 150 63 INTEGER, PARAMETER :: nb_point_max = 2000 64 INTEGER, PARAMETER :: nb_type_class = 10 65 INTEGER, PARAMETER :: nb_3d_vars = 3 66 INTEGER, PARAMETER :: nb_2d_vars = 2 67 INTEGER :: nb_sec 68 69 TYPE POINT_SECTION 70 INTEGER :: I,J 71 END TYPE POINT_SECTION 72 73 TYPE COORD_SECTION 74 REAL(wp) :: lon,lat 75 END TYPE COORD_SECTION 76 77 TYPE SECTION 78 CHARACTER(len=60) :: name ! name of the sec 79 LOGICAL :: llstrpond ! true if you want the computation of salt and 91 80 ! heat transports 92 LOGICAL :: ll_ice_section ! ice surface and ice volume computation 93 LOGICAL :: ll_date_line ! = T if the section crosses the date-line 94 TYPE(COORD_SECTION), DIMENSION(2) :: coordSec ! longitude and latitude of the extremities of the sec 95 INTEGER :: nb_class ! number of boundaries for density classes 96 INTEGER, DIMENSION(nb_point_max) :: direction ! vector direction of the point in the section 97 CHARACTER(len=40),DIMENSION(nb_class_max) :: classname ! characteristics of the class 98 REAL(wp), DIMENSION(nb_class_max) :: zsigi ,&! in-situ density classes (99 if you don't want) 99 zsigp ,&! potential density classes (99 if you don't want) 100 zsal ,&! salinity classes (99 if you don't want) 101 ztem ,&! temperature classes(99 if you don't want) 102 zlay ! level classes (99 if you don't want) 103 REAL(wp), DIMENSION(nb_type_class,nb_class_max) :: transport ! transport output 104 REAL(wp) :: slopeSection ! slope of the section 105 INTEGER :: nb_point ! number of points in the section 106 TYPE(POINT_SECTION),DIMENSION(nb_point_max) :: listPoint ! list of points in the sections 107 END TYPE SECTION 108 109 TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections 110 111 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: transports_3d 112 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: transports_2d 113 81 LOGICAL :: ll_ice_section ! ice surface and ice volume computation 82 LOGICAL :: ll_date_line ! = T if the section crosses the date-line 83 TYPE(COORD_SECTION), DIMENSION(2) :: coordSec ! longitude and latitude of the extremities of the sec 84 INTEGER :: nb_class ! number of boundaries for density classes 85 INTEGER, DIMENSION(nb_point_max) :: direction ! vector direction of the point in the section 86 CHARACTER(len=40),DIMENSION(nb_class_max) :: classname ! characteristics of the class 87 REAL(wp), DIMENSION(nb_class_max) :: zsigi ,&! in-situ density classes (99 if you don't want) 88 zsigp ,&! potential density classes (99 if you don't want) 89 zsal ,&! salinity classes (99 if you don't want) 90 ztem ,&! temperature classes(99 if you don't want) 91 zlay ! level classes (99 if you don't want) 92 REAL(wp), DIMENSION(nb_type_class,nb_class_max) :: transport ! transport output 93 REAL(wp) :: slopeSection ! slope of the section 94 INTEGER :: nb_point ! number of points in the section 95 TYPE(POINT_SECTION),DIMENSION(nb_point_max) :: listPoint ! list of points in the sections 96 END TYPE SECTION 97 98 TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections 99 100 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: transports_3d 101 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: transports_2d 102 103 !!---------------------------------------------------------------------- 104 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 114 105 !! $Id$ 106 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 107 !!---------------------------------------------------------------------- 115 108 CONTAINS 116 117 109 118 110 INTEGER FUNCTION diadct_alloc() … … 130 122 131 123 END FUNCTION diadct_alloc 124 132 125 133 126 SUBROUTINE dia_dct_init … … 208 201 !! Reinitialise all relevant arrays to zero 209 202 !!--------------------------------------------------------------------- 210 !! * Arguments 211 INTEGER,INTENT(IN) ::kt 212 213 !! * Local variables 203 INTEGER,INTENT(in) ::kt 204 ! 214 205 INTEGER :: jsec, &! loop on sections 215 206 itotal ! nb_sec_max*nb_type_class*nb_class_max … … 220 211 REAL(wp), POINTER, DIMENSION(:) :: zwork ! " 221 212 REAL(wp), POINTER, DIMENSION(:,:,:):: zsum ! " 222 223 213 !!--------------------------------------------------------------------- 214 ! 224 215 IF( nn_timing == 1 ) CALL timing_start('dia_dct') 225 216 … … 619 610 zumid_ice, zvmid_ice, &!U/V ice velocity 620 611 zTnorm !transport of velocity through one cell's sides 621 REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, z fsdep !temperature/salinity/potential density/ssh/depth at u/v point612 REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep !temperature/salinity/potential density/ssh/depth at u/v point 622 613 623 614 TYPE(POINT_SECTION) :: k 624 !!--------------------------------------------------------625 626 IF( ld_debug )WRITE(numout,*)' Compute transport'627 628 !---------------------------!629 ! COMPUTE TRANSPORT !630 !---------------------------!631 IF(sec%nb_point .NE. 0)THEN632 633 !----------------------------------------------------------------------------------------------------634 !Compute sign for velocities:635 !636 !convention:637 ! non horizontal section: direction + is toward left hand of section638 ! horizontal section: direction + is toward north of section639 !640 !641 ! slopeSection < 0 slopeSection > 0 slopeSection=inf slopeSection=0642 ! ---------------- ----------------- --------------- --------------643 !644 ! isgnv=1 direction +645 ! ______ _____ ______646 ! | //| | | direction +647 ! | isgnu=1 // | |isgnu=1 |isgnu=1 /|\648 ! |_______ // ______| \\ | ---\ |649 ! | | isgnv=-1 \\ | | ---/ direction + ____________650 ! | | __\\| |651 ! | | direction + | isgnv=1652 !653 !----------------------------------------------------------------------------------------------------654 isgnu = 1655 IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1656 ELSE ; isgnv = 1657 ENDIF658 IF( sec%slopeSection .GE. 9999. ) isgnv = 1659 660 IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv661 662 !--------------------------------------!663 ! LOOP ON THE SEGMENT BETWEEN 2 NODES !664 !--------------------------------------!665 DO jseg=1,MAX(sec%nb_point-1,0)615 !!-------------------------------------------------------- 616 ! 617 IF( ld_debug )WRITE(numout,*)' Compute transport' 618 619 !---------------------------! 620 ! COMPUTE TRANSPORT ! 621 !---------------------------! 622 IF(sec%nb_point .NE. 0)THEN 623 624 !---------------------------------------------------------------------------------------------------- 625 !Compute sign for velocities: 626 ! 627 !convention: 628 ! non horizontal section: direction + is toward left hand of section 629 ! horizontal section: direction + is toward north of section 630 ! 631 ! 632 ! slopeSection < 0 slopeSection > 0 slopeSection=inf slopeSection=0 633 ! ---------------- ----------------- --------------- -------------- 634 ! 635 ! isgnv=1 direction + 636 ! ______ _____ ______ 637 ! | //| | | direction + 638 ! | isgnu=1 // | |isgnu=1 |isgnu=1 /|\ 639 ! |_______ // ______| \\ | ---\ | 640 ! | | isgnv=-1 \\ | | ---/ direction + ____________ 641 ! | | __\\| | 642 ! | | direction + | isgnv=1 643 ! 644 !---------------------------------------------------------------------------------------------------- 645 isgnu = 1 646 IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1 647 ELSE ; isgnv = 1 648 ENDIF 649 IF( sec%slopeSection .GE. 9999. ) isgnv = 1 650 651 IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv 652 653 !--------------------------------------! 654 ! LOOP ON THE SEGMENT BETWEEN 2 NODES ! 655 !--------------------------------------! 656 DO jseg=1,MAX(sec%nb_point-1,0) 666 657 667 !------------------------------------------------------------------------------------------- 668 ! Select the appropriate coordinate for computing the velocity of the segment 669 ! 670 ! CASE(0) Case (2) 671 ! ------- -------- 672 ! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j) 673 ! F(i,j)----------V(i+1,j)-------F(i+1,j) | 674 ! | 675 ! | 676 ! | 677 ! Case (3) U(i,j) 678 ! -------- | 679 ! | 680 ! listPoint(jseg+1) F(i,j+1) | 681 ! | | 682 ! | | 683 ! | listPoint(jseg+1) F(i,j-1) 684 ! | 685 ! | 686 ! U(i,j+1) 687 ! | Case(1) 688 ! | ------ 689 ! | 690 ! | listPoint(jseg+1) listPoint(jseg) 691 ! | F(i-1,j)-----------V(i,j) -------f(jseg) 692 ! listPoint(jseg) F(i,j) 693 ! 694 !------------------------------------------------------------------------------------------- 695 696 SELECT CASE( sec%direction(jseg) ) 697 CASE(0) ; k = sec%listPoint(jseg) 698 CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 699 CASE(2) ; k = sec%listPoint(jseg) 700 CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 701 END SELECT 702 703 !---------------------------| 704 ! LOOP ON THE LEVEL | 705 !---------------------------| 706 !Sum of the transport on the vertical 707 DO jk=1,mbathy(k%I,k%J) 708 709 ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 710 SELECT CASE( sec%direction(jseg) ) 711 CASE(0,1) 712 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 713 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 714 zrhop = interp(k%I,k%J,jk,'V',rhop) 715 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 716 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 717 CASE(2,3) 718 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 719 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 720 zrhop = interp(k%I,k%J,jk,'U',rhop) 721 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 722 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 723 END SELECT 724 725 zfsdep= fsdept(k%I,k%J,jk) 658 !------------------------------------------------------------------------------------------- 659 ! Select the appropriate coordinate for computing the velocity of the segment 660 ! 661 ! CASE(0) Case (2) 662 ! ------- -------- 663 ! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j) 664 ! F(i,j)----------V(i+1,j)-------F(i+1,j) | 665 ! | 666 ! | 667 ! | 668 ! Case (3) U(i,j) 669 ! -------- | 670 ! | 671 ! listPoint(jseg+1) F(i,j+1) | 672 ! | | 673 ! | | 674 ! | listPoint(jseg+1) F(i,j-1) 675 ! | 676 ! | 677 ! U(i,j+1) 678 ! | Case(1) 679 ! | ------ 680 ! | 681 ! | listPoint(jseg+1) listPoint(jseg) 682 ! | F(i-1,j)-----------V(i,j) -------f(jseg) 683 ! listPoint(jseg) F(i,j) 684 ! 685 !------------------------------------------------------------------------------------------- 686 687 SELECT CASE( sec%direction(jseg) ) 688 CASE(0) ; k = sec%listPoint(jseg) 689 CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 690 CASE(2) ; k = sec%listPoint(jseg) 691 CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 692 END SELECT 693 694 !---------------------------| 695 ! LOOP ON THE LEVEL | 696 !---------------------------| 697 DO jk = 1, mbathy(k%I,k%J) !Sum of the transport on the vertical 698 ! ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 699 SELECT CASE( sec%direction(jseg) ) 700 CASE(0,1) 701 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 702 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 703 zrhop = interp(k%I,k%J,jk,'V',rhop) 704 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 705 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 706 CASE(2,3) 707 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 708 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 709 zrhop = interp(k%I,k%J,jk,'U',rhop) 710 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 711 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 712 END SELECT 713 ! 714 zdep= gdept_n(k%I,k%J,jk) 726 715 727 !compute velocity with the correct direction728 SELECT CASE( sec%direction(jseg) )729 CASE(0,1)730 zumid=0.731 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)732 CASE(2,3)733 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)734 zvmid=0.735 END SELECT736 737 !zTnorm=transport through one cell;738 !velocity* cell's length * cell's thickness739 zTnorm=zumid*e2u(k%I,k%J)* fse3u(k%I,k%J,jk)+ &740 zvmid*e1v(k%I,k%J)* fse3v(k%I,k%J,jk) 741 742 #if ! defined key_vvl 743 !add transport due to free surface744 IF( jk==1 )THEN745 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &746 zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)747 ENDIF748 #endif 716 SELECT CASE( sec%direction(jseg) ) !compute velocity with the correct direction 717 CASE(0,1) 718 zumid=0._wp 719 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 720 CASE(2,3) 721 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 722 zvmid=0._wp 723 END SELECT 724 725 !zTnorm=transport through one cell; 726 !velocity* cell's length * cell's thickness 727 zTnorm = zumid*e2u(k%I,k%J) * e3u_n(k%I,k%J,jk) & 728 & + zvmid*e1v(k%I,k%J) * e3v_n(k%I,k%J,jk) 729 730 !!gm THIS is WRONG no transport due to ssh in linear free surface case !!!!! 731 IF( ln_linssh ) THEN !add transport due to free surface 732 IF( jk==1 ) THEN 733 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) & 734 & + zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 735 ENDIF 736 ENDIF 737 !!gm end 749 738 !COMPUTE TRANSPORT 750 739 751 740 transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 752 741 753 IF 742 IF( sec%llstrpond ) THEN 754 743 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * ztn * zrhop * rcp 755 744 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zsn * zrhop * 0.001 756 745 ENDIF 757 746 758 END DO !end of loop on the level747 END DO !end of loop on the level 759 748 760 749 #if defined key_lim2 || defined key_lim3 … … 797 786 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* & 798 787 a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) 799 END DO788 END DO 800 789 #endif 801 790 … … 803 792 #endif 804 793 805 END DO !end of loop on the segment794 END DO !end of loop on the segment 806 795 807 796 ENDIF !end of sec%nb_point =0 case 808 797 ! 809 798 END SUBROUTINE transport 810 799 800 811 801 SUBROUTINE dia_dct_sum(sec,jsec) 812 802 !!------------------------------------------------------------- … … 828 818 !! 829 819 !!------------------------------------------------------------- 830 !! * arguments831 820 TYPE(SECTION),INTENT(INOUT) :: sec 832 821 INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section … … 834 823 TYPE(POINT_SECTION) :: k 835 824 INTEGER :: jk,jseg,jclass ! dummy variables for looping on level/segment/classes 836 REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, z fsdep ! temperature/salinity/ssh/potential density /depth at u/v point825 REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/ssh/potential density /depth at u/v point 837 826 !!------------------------------------------------------------- 838 827 … … 903 892 END SELECT 904 893 905 z fsdep= fsdept(k%I,k%J,jk)894 zdep= gdept_n(k%I,k%J,jk) 906 895 907 896 !------------------------------- … … 932 921 ( sec%ztem(jclass) .EQ.99.)) .AND. & 933 922 934 ((( z fsdep .GE. sec%zlay(jclass)) .AND. &935 ( z fsdep .LE. sec%zlay(jclass+1))) .OR. &923 ((( zdep .GE. sec%zlay(jclass)) .AND. & 924 ( zdep .LE. sec%zlay(jclass+1))) .OR. & 936 925 ( sec%zlay(jclass) .EQ. 99. )) & 937 926 )) THEN … … 991 980 #endif 992 981 993 END DO !end of loop on the segment982 END DO !end of loop on the segment 994 983 995 984 ELSE !if sec%nb_point =0 … … 1000 989 1001 990 END SUBROUTINE dia_dct_sum 1002 991 992 1003 993 SUBROUTINE dia_dct_wri(kt,ksec,sec) 1004 994 !!------------------------------------------------------------- … … 1138 1128 sec%transport(9,1),sec%transport(10,1), & 1139 1129 sec%transport(9,1)+sec%transport(10,1) 1140 ENDIF1130 ENDIF 1141 1131 1142 118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4) 1143 119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6) 1144 1145 CALL wrk_dealloc(nb_type_class , zsumclasses ) 1146 END SUBROUTINE dia_dct_wri 1147 1148 FUNCTION interp(ki, kj, kk, cd_point, ptab) 1132 118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4) 1133 119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6) 1134 1135 CALL wrk_dealloc(nb_type_class , zsumclasses ) 1136 ! 1137 END SUBROUTINE dia_dct_wri 1138 1139 1140 FUNCTION interp(ki, kj, kk, cd_point, ptab) 1149 1141 !!---------------------------------------------------------------------- 1150 1142 !! … … 1214 1206 !*local declations 1215 1207 INTEGER :: ii1, ij1, ii2, ij2 ! local integer 1216 REAL(wp):: ze3t, z fse3, zwgt1, zwgt2, zbis, zdepu ! local real1208 REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu ! local real 1217 1209 REAL(wp):: zet1, zet2 ! weight for interpolation 1218 1210 REAL(wp):: zdep1,zdep2 ! differences of depth … … 1241 1233 IF( ln_sco )THEN ! s-coordinate case 1242 1234 1243 zdepu = ( fsdept(ii1,ij1,kk) + fsdept(ii2,ij2,kk) ) /21244 zdep1 = fsdept(ii1,ij1,kk) - zdepu1245 zdep2 = fsdept(ii2,ij2,kk) - zdepu1235 zdepu = ( gdept_n(ii1,ij1,kk) + gdept_n(ii2,ij2,kk) ) * 0.5_wp 1236 zdep1 = gdept_n(ii1,ij1,kk) - zdepu 1237 zdep2 = gdept_n(ii2,ij2,kk) - zdepu 1246 1238 1247 1239 ! weights … … 1255 1247 ELSE ! full step or partial step case 1256 1248 1257 #if defined key_vvl 1258 1259 ze3t = fse3t_n(ii2,ij2,kk) - fse3t_n(ii1,ij1,kk) 1260 zwgt1 = ( fse3w_n(ii2,ij2,kk) - fse3w_n(ii1,ij1,kk) ) / fse3w_n(ii2,ij2,kk) 1261 zwgt2 = ( fse3w_n(ii1,ij1,kk) - fse3w_n(ii2,ij2,kk) ) / fse3w_n(ii1,ij1,kk) 1262 1263 #else 1264 1265 ze3t = fse3t(ii2,ij2,kk) - fse3t(ii1,ij1,kk) 1266 zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk) 1267 zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk) 1268 1269 #endif 1249 ze3t = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk) 1250 zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk) 1251 zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk) 1270 1252 1271 1253 IF(kk .NE. 1)THEN … … 1288 1270 1289 1271 ENDIF 1290 1291 1292 END FUNCTION interp 1272 ! 1273 END FUNCTION interp 1293 1274 1294 1275 #else … … 1311 1292 #endif 1312 1293 1294 !!====================================================================== 1313 1295 END MODULE diadct -
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90
r5836 r6140 33 33 34 34 !! * Substitutions 35 # include "domzgr_substitute.h90"36 35 # include "vectopt_loop_substitute.h90" 37 36 !!---------------------------------------------------------------------- … … 40 39 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 41 40 !!---------------------------------------------------------------------- 42 43 41 CONTAINS 44 42 … … 80 78 DO jj = 2, jpjm1 81 79 DO ji = fs_2, fs_jpim1 ! vector opt. 82 zwei = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)80 zwei = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 83 81 a_salb = a_salb + ( tsb(ji,jj,jk,jp_sal) - zsm0 ) * zwei 84 82 END DO … … 106 104 DO jj = 2, jpjm1 107 105 DO ji = fs_2, fs_jpim1 ! vector opt. 108 zwei = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)106 zwei = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 109 107 a_saln = a_saln + ( tsn(ji,jj,jk,jp_sal) - zsm0 ) * zwei 110 108 zvol = zvol + zwei … … 116 114 117 115 ! Conversion in m3 118 a_fwf = a_fwf * rdt tra(1)* 1.e-3116 a_fwf = a_fwf * rdt * 1.e-3 119 117 120 118 ! fwf correction to bring back the mean ssh to zero … … 186 184 zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 187 185 zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 188 zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)186 zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj) 189 187 190 188 IF( un(ji,jj,jk) > 0.e0 ) THEN … … 238 236 zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 239 237 zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 240 zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)238 zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj) 241 239 242 240 IF( un(ji,jj,jk) > 0.e0 ) THEN … … 290 288 zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 291 289 zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 292 zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)290 zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj) 293 291 294 292 IF( un(ji,jj,jk) > 0.e0 ) THEN … … 342 340 zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 343 341 zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 344 zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)342 zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj) 345 343 346 344 IF( un(ji,jj,jk) > 0.e0 ) THEN … … 404 402 WRITE(inum,*) 405 403 WRITE(inum,*) 'Net freshwater budget ' 406 WRITE(inum,9010) ' fwf = ',a_fwf, ' m3 =', a_fwf /(FLOAT(nitend-nit000+1)*rdt tra(1)) * 1.e-6,' Sv'404 WRITE(inum,9010) ' fwf = ',a_fwf, ' m3 =', a_fwf /(FLOAT(nitend-nit000+1)*rdt) * 1.e-6,' Sv' 407 405 WRITE(inum,*) 408 406 WRITE(inum,9010) ' zarea =',zarea … … 460 458 ENDIF 461 459 462 IF( nn_timing == 1 ) CALL timing_st art('dia_fwb')460 IF( nn_timing == 1 ) CALL timing_stop('dia_fwb') 463 461 464 462 9005 FORMAT(1X,A,ES24.16) -
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
r5930 r6140 21 21 USE ioipsl ! NetCDF IPSL library 22 22 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 23 USE diadimg ! To write dimg24 23 USE timing ! preformance summary 25 24 USE wrk_nemo ! working arrays … … 135 134 DO jk=1,nb_ana 136 135 DO ji=1,jpmax_harmo 137 IF (TRIM(tname(jk)) .eq.Wave(ji)%cname_tide) THEN136 IF (TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN 138 137 name(jk) = ji 139 138 EXIT … … 194 193 DO ji = 1,jpi 195 194 ! Elevation 196 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)* tmask_i(ji,jj)197 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)* umask_i(ji,jj)198 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)* vmask_i(ji,jj)195 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj) 196 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 197 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 199 198 END DO 200 199 END DO … … 324 323 X1= ana_amp(ji,jj,jh,1) 325 324 X2=-ana_amp(ji,jj,jh,2) 326 out_u(ji,jj, jh) = X1 * umask_i(ji,jj)327 out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)325 out_u(ji,jj, jh) = X1 * ssumask(ji,jj) 326 out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj) 328 327 ENDDO 329 328 ENDDO … … 358 357 X1=ana_amp(ji,jj,jh,1) 359 358 X2=-ana_amp(ji,jj,jh,2) 360 out_v(ji,jj, jh)=X1 * vmask_i(ji,jj)361 out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)359 out_v(ji,jj, jh)=X1 * ssvmask(ji,jj) 360 out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj) 362 361 END DO 363 362 END DO … … 384 383 !!---------------------------------------------------------------------- 385 384 386 #if defined key_dimgout387 cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'388 cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'389 cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'390 #endif391 392 385 IF(lwp) WRITE(numout,*) ' ' 393 386 IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results' 394 #if defined key_dimgout395 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ Output files: ', TRIM(cdfile_name_T)396 IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_U)397 IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_V)398 #endif399 387 IF(lwp) WRITE(numout,*) ' ' 400 388 … … 402 390 !///////////// 403 391 ! 404 #if defined key_dimgout405 cltext='Elevation amplitude and phase'406 CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')407 #else408 392 DO jh = 1, nb_ana 409 393 CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) ) 410 394 CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) ) 411 395 END DO 412 #endif413 396 414 397 ! B) ubar 415 398 !///////// 416 399 ! 417 #if defined key_dimgout418 cltext='ubar amplitude and phase'419 CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')420 #else421 400 DO jh = 1, nb_ana 422 401 CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) ) 423 402 CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) ) 424 403 END DO 425 #endif426 404 427 405 ! C) vbar 428 406 !///////// 429 407 ! 430 #if defined key_dimgout431 cltext='vbar amplitude and phase'432 CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')433 #else434 408 DO jh = 1, nb_ana 435 409 CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh ) ) 436 410 CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) ) 437 411 END DO 438 #endif439 412 ! 440 413 END SUBROUTINE dia_wri_harm … … 488 461 DO jj_sd = ji_sd, ninco 489 462 zval2 = ABS(ztmp3(ji_sd,jj_sd)) 490 IF( zval2 .GE.zval1 )THEN463 IF( zval2 >= zval1 )THEN 491 464 ipivot(ji_sd) = jj_sd 492 465 zval1 = zval2 -
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
r5643 r6140 46 46 REAL(wp) :: frc_wn_t, frc_wn_s ! global forcing trends 47 47 ! 48 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf , ssh_ini ! 48 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf 49 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf_ini , ssh_ini ! 49 50 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: ssh_hc_loc_ini, ssh_sc_loc_ini ! 50 51 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: hc_loc_ini, sc_loc_ini, e3t_ini ! 51 52 52 53 !! * Substitutions 53 # include "domzgr_substitute.h90"54 54 # include "vectopt_loop_substitute.h90" 55 55 !!---------------------------------------------------------------------- … … 100 100 IF( ln_rnf_sal) z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) ) 101 101 ! Add ice shelf heat & salt input 102 IF( nn_isf .GE. 1 ) THEN 103 z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) ) 104 z_frc_trd_s = z_frc_trd_s + glob_sum( risf_tsc(:,:,jp_sal) * surf(:,:) ) 105 ENDIF 106 102 IF( ln_isf ) z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) ) 107 103 ! Add penetrative solar radiation 108 104 IF( ln_traqsr ) z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr (:,:) * surf(:,:) ) … … 110 106 IF( ln_trabbc ) z_frc_trd_t = z_frc_trd_t + glob_sum( qgh_trd0(:,:) * surf(:,:) ) 111 107 ! 112 IF( .NOT. lk_vvl) THEN113 IF 108 IF( ln_linssh ) THEN 109 IF( ln_isfcav ) THEN 114 110 DO ji=1,jpi 115 111 DO jj=1,jpj 116 112 z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem) 117 113 z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal) 118 END DO119 END DO114 END DO 115 END DO 120 116 ELSE 121 117 z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) … … 130 126 frc_s = frc_s + z_frc_trd_s * rdt 131 127 ! ! Advection flux through fixed surface (z=0) 132 IF( .NOT. lk_vvl) THEN128 IF( ln_linssh ) THEN 133 129 frc_wn_t = frc_wn_t + z_wn_trd_t * rdt 134 130 frc_wn_s = frc_wn_s + z_wn_trd_s * rdt … … 138 134 ! 2 - Content variations ! 139 135 ! ------------------------ ! 136 ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl) 140 137 zdiff_v2 = 0._wp 141 138 zdiff_hc = 0._wp … … 143 140 144 141 ! volume variation (calculated with ssh) 145 zdiff_v1 = glob_sum ( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:)) )142 zdiff_v1 = glob_sum_full( surf(:,:) * sshn(:,:) - surf_ini(:,:) * ssh_ini(:,:) ) 146 143 147 144 ! heat & salt content variation (associated with ssh) 148 IF( .NOT. lk_vvl) THEN149 IF 145 IF( ln_linssh ) THEN 146 IF( ln_isfcav ) THEN 150 147 DO ji = 1, jpi 151 148 DO jj = 1, jpj … … 158 155 z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) 159 156 END IF 160 z_ssh_hc = glob_sum ( z2d0 )161 z_ssh_sc = glob_sum ( z2d1 )157 z_ssh_hc = glob_sum_full( z2d0 ) 158 z_ssh_sc = glob_sum_full( z2d1 ) 162 159 ENDIF 163 160 164 161 DO jk = 1, jpkm1 165 162 ! volume variation (calculated with scale factors) 166 zdiff_v2 = zdiff_v2 + glob_sum ( surf(:,:) * tmask(:,:,jk)&167 & * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk)) )163 zdiff_v2 = zdiff_v2 + glob_sum_full( surf(:,:) * tmask(:,:,jk) & 164 & * e3t_n(:,:,jk) - surf_ini(:,:) * e3t_ini(:,:,jk) ) 168 165 ! heat content variation 169 zdiff_hc = zdiff_hc + glob_sum ( surf(:,:) * tmask(:,:,jk)&170 & * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) )166 zdiff_hc = zdiff_hc + glob_sum_full( surf(:,:) * tmask(:,:,jk) & 167 & * e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - surf_ini(:,:) * hc_loc_ini(:,:,jk) ) 171 168 ! salt content variation 172 zdiff_sc = zdiff_sc + glob_sum ( surf(:,:) * tmask(:,:,jk)&173 & * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk)) )169 zdiff_sc = zdiff_sc + glob_sum_full( surf (:,:) * tmask(:,:,jk) & 170 * e3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - surf_ini(:,:) * sc_loc_ini(:,:,jk) ) 174 171 ENDDO 175 172 176 173 ! Substract forcing from heat content, salt content and volume variations 177 174 zdiff_v1 = zdiff_v1 - frc_v 178 IF( lk_vvl) zdiff_v2 = zdiff_v2 - frc_v175 IF( .NOT.ln_linssh ) zdiff_v2 = zdiff_v2 - frc_v 179 176 zdiff_hc = zdiff_hc - frc_t 180 177 zdiff_sc = zdiff_sc - frc_s 181 IF( .NOT. lk_vvl) THEN178 IF( ln_linssh ) THEN 182 179 zdiff_hc1 = zdiff_hc + z_ssh_hc 183 180 zdiff_sc1 = zdiff_sc + z_ssh_sc … … 191 188 zvol_tot = 0._wp ! total ocean volume (calculated with scale factors) 192 189 DO jk = 1, jpkm1 193 zvol_tot = zvol_tot + glob_sum ( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )190 zvol_tot = zvol_tot + glob_sum_full( surf(:,:) * tmask(:,:,jk) * e3t_n(:,:,jk) ) 194 191 END DO 195 192 196 193 !!gm to be added ? 197 ! IF( .NOT. lk_vvl) THEN ! fixed volume, add the ssh contribution194 ! IF( ln_linssh ) THEN ! fixed volume, add the ssh contribution 198 195 ! zvol_tot = zvol_tot + glob_sum( surf(:,:) * sshn(:,:) ) 199 196 ! ENDIF 200 197 !!gm end 201 198 202 IF( lk_vvl ) THEN 199 IF( ln_linssh ) THEN 200 CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot) ! Heat content variation (C) 201 CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot) ! Salt content variation (psu) 202 CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp ) ! Heat content variation (1.e20 J) 203 CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9 ) ! Salt content variation (psu*km3) 204 CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 ) ! volume ssh variation (km3) 205 CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 ) ! vol - surface forcing (km3) 206 CALL iom_put( 'bgfrctem' , frc_t / zvol_tot ) ! hc - surface forcing (C) 207 CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot ) ! sc - surface forcing (psu) 208 CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot ) ! hc - error due to free surface (C) 209 CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot ) ! sc - error due to free surface (psu) 210 ELSE 203 211 CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot ) ! Temperature variation (C) 204 212 CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot ) ! Salinity variation (psu) … … 210 218 CALL iom_put( 'bgfrctem' , frc_t / zvol_tot ) ! hc - surface forcing (C) 211 219 CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot ) ! sc - surface forcing (psu) 212 ELSE213 CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot) ! Heat content variation (C)214 CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot) ! Salt content variation (psu)215 CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp ) ! Heat content variation (1.e20 J)216 CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9 ) ! Salt content variation (psu*km3)217 CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 ) ! volume ssh variation (km3)218 CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 ) ! vol - surface forcing (km3)219 CALL iom_put( 'bgfrctem' , frc_t / zvol_tot ) ! hc - surface forcing (C)220 CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot ) ! sc - surface forcing (psu)221 CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot ) ! hc - error due to free surface (C)222 CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot ) ! sc - error due to free surface (psu)223 220 ENDIF 224 221 ! 225 222 IF( lrst_oce ) CALL dia_hsb_rst( kt, 'WRITE' ) 226 223 ! 227 224 CALL wrk_dealloc( jpi,jpj, z2d0, z2d1 ) 228 225 ! 229 226 IF( nn_timing == 1 ) CALL timing_stop('dia_hsb') 230 227 ! … … 257 254 CALL iom_get( numror, 'frc_t', frc_t ) 258 255 CALL iom_get( numror, 'frc_s', frc_s ) 259 IF( .NOT. lk_vvl) THEN256 IF( ln_linssh ) THEN 260 257 CALL iom_get( numror, 'frc_wn_t', frc_wn_t ) 261 258 CALL iom_get( numror, 'frc_wn_s', frc_wn_s ) 262 259 ENDIF 260 CALL iom_get( numror, jpdom_autoglo, 'surf_ini', surf_ini ) ! ice sheet coupling 263 261 CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini ) 264 262 CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini ) 265 263 CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini ) 266 264 CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini ) 267 IF( .NOT. lk_vvl) THEN265 IF( ln_linssh ) THEN 268 266 CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini ) 269 267 CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini ) … … 273 271 IF(lwp) WRITE(numout,*) ' dia_hsb at initial state ' 274 272 IF(lwp) WRITE(numout,*) '~~~~~~~' 275 ssh_ini(:,:) = sshn(:,:) ! initial ssh 273 surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:) ! initial ocean surface 274 ssh_ini(:,:) = sshn(:,:) ! initial ssh 276 275 DO jk = 1, jpk 277 e3t_ini (:,:,jk) = fse3t_n(:,:,jk) ! initial vertical scale factors 278 hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) ! initial heat content 279 sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) ! initial salt content 276 ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance). 277 e3t_ini (:,:,jk) = e3t_n(:,:,jk) * tmask(:,:,jk) ! initial vertical scale factors 278 hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * e3t_n(:,:,jk) * tmask(:,:,jk) ! initial heat content 279 sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * e3t_n(:,:,jk) * tmask(:,:,jk) ! initial salt content 280 280 END DO 281 281 frc_v = 0._wp ! volume trend due to forcing 282 282 frc_t = 0._wp ! heat content - - - - 283 283 frc_s = 0._wp ! salt content - - - - 284 IF( .NOT. lk_vvl) THEN284 IF( ln_linssh ) THEN 285 285 IF ( ln_isfcav ) THEN 286 286 DO ji=1,jpi … … 308 308 CALL iom_rstput( kt, nitrst, numrow, 'frc_t' , frc_t ) 309 309 CALL iom_rstput( kt, nitrst, numrow, 'frc_s' , frc_s ) 310 IF( .NOT. lk_vvl) THEN310 IF( ln_linssh ) THEN 311 311 CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t ) 312 312 CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s ) 313 313 ENDIF 314 CALL iom_rstput( kt, nitrst, numrow, 'surf_ini', surf_ini ) ! ice sheet coupling 314 315 CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini ) 315 316 CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini ) 316 317 CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini ) 317 318 CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini ) 318 IF( .NOT. lk_vvl) THEN319 IF( ln_linssh ) THEN 319 320 CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini ) 320 321 CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini ) … … 379 380 ! 1 - Allocate memory ! 380 381 ! ------------------- ! 381 ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), &382 & e3t_ini(jpi,jpj,jpk), surf(jpi,jpj), ssh_ini(jpi,jpj), STAT=ierror )382 ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), & 383 & e3t_ini(jpi,jpj,jpk), surf(jpi,jpj), ssh_ini(jpi,jpj), STAT=ierror ) 383 384 IF( ierror > 0 ) THEN 384 385 CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) ; RETURN 385 386 ENDIF 386 387 387 IF( .NOT. lk_vvl )ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )388 IF( ln_linssh ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror ) 388 389 IF( ierror > 0 ) THEN 389 390 CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) ; RETURN -
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90
r5836 r6140 20 20 USE dom_oce ! ocean space and time domain 21 21 USE phycst ! physical constants 22 ! 22 23 USE in_out_manager ! I/O manager 23 24 USE lib_mpp ! MPP library … … 31 32 PUBLIC dia_hth_alloc ! routine called by nemogcm.F90 32 33 33 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .TRUE. !: thermocline-20d depths flag 34 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .TRUE. !: thermocline-20d depths flag 35 34 36 ! note: following variables should move to local variables once iom_put is always used 35 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hth !: depth of the max vertical temperature gradient [m] … … 38 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htc3 !: heat content of first 300 m [W] 39 41 40 !! * Substitutions41 # include "domzgr_substitute.h90"42 42 !!---------------------------------------------------------------------- 43 43 !! NEMO/OPA 4.0 , NEMO Consortium (2011) … … 52 52 !!--------------------------------------------------------------------- 53 53 ! 54 ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc)54 ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc ) 55 55 ! 56 56 IF( lk_mpp ) CALL mpp_sum ( dia_hth_alloc ) … … 108 108 IF( kt == nit000 ) THEN 109 109 ! ! allocate dia_hth array 110 IF( dia_hth_alloc() /= 0 ) CALL ctl_stop( 'STOP', ' lim_sbc_init: unable to allocate standard arrays' )111 112 IF(. not. ALLOCATED(ik20))THEN110 IF( dia_hth_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' ) 111 112 IF(.NOT. ALLOCATED(ik20) ) THEN 113 113 ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), & 114 114 & zabs2(jpi,jpj), & … … 187 187 DO ji = 1, jpi 188 188 ! 189 zzdep = fsdepw(ji,jj,jk)189 zzdep = gdepw_n(ji,jj,jk) 190 190 zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) 191 191 zzdep = zzdep * tmask(ji,jj,1) … … 223 223 DO ji = 1, jpi 224 224 ! 225 zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1)225 zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 226 226 ! 227 227 zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem) ! - delta T(10m) … … 270 270 DO ji = 1, jpi 271 271 ! 272 zzdep = fsdepw(ji,jj,mbkt(ji,jj)+1) ! depth of the oean bottom272 zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1) ! depth of the oean bottom 273 273 ! 274 274 iid = ik20(ji,jj) 275 275 IF( iid /= 1 ) THEN 276 zztmp = fsdept(ji,jj,iid ) & ! linear interpolation277 & + ( fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) ) &276 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 277 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 278 278 & * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 279 279 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) … … 285 285 iid = ik28(ji,jj) 286 286 IF( iid /= 1 ) THEN 287 zztmp = fsdept(ji,jj,iid ) & ! linear interpolation288 & + ( fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) ) &287 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 288 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 289 289 & * ( 28.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 290 290 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) … … 311 311 END DO 312 312 ! surface boundary condition 313 IF( l k_vvl ) THEN ; zthick(:,:) = 0._wp ; htc3(:,:) = 0._wp314 ELSE ; zthick(:,:) = sshn(:,:) ; htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)313 IF( ln_linssh ) THEN ; zthick(:,:) = sshn(:,:) ; htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 314 ELSE ; zthick(:,:) = 0._wp ; htc3(:,:) = 0._wp 315 315 ENDIF 316 316 ! integration down to ilevel 317 317 DO jk = 1, ilevel 318 zthick(:,:) = zthick(:,:) + fse3t(:,:,jk)319 htc3 (:,:) = htc3 (:,:) + fse3t(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)318 zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 319 htc3 (:,:) = htc3 (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 320 320 END DO 321 321 ! deepest layer … … 323 323 DO jj = 1, jpj 324 324 DO ji = 1, jpi 325 htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) )&326 325 htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) & 326 & * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 327 327 END DO 328 328 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/dianam.F90
r2528 r6140 72 72 73 73 IF( llfsec .OR. kfreq < 0 ) THEN ; inbsec = kfreq ! output frequency already in seconds 74 ELSE ; inbsec = kfreq * NINT( rdt tra(1)) ! from time-step to seconds74 ELSE ; inbsec = kfreq * NINT( rdt ) ! from time-step to seconds 75 75 ENDIF 76 76 iddss = NINT( rday ) ! number of seconds in 1 day … … 116 116 ! date of the beginning and the end of the run 117 117 118 zdrun = rdt tra(1)/ rday * REAL( nitend - nit000, wp ) ! length of the run in days119 zjul = fjulday - rdt tra(1)/ rday118 zdrun = rdt / rday * REAL( nitend - nit000, wp ) ! length of the run in days 119 zjul = fjulday - rdt / rday 120 120 CALL ju2ymds( zjul , iyear1, imonth1, iday1, zsec1 ) ! year/month/day of the beginning of run 121 121 CALL ju2ymds( zjul + zdrun, iyear2, imonth2, iday2, zsec2 ) ! year/month/day of the end of run -
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
r5147 r6140 59 59 REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:) :: p_fval2d 60 60 61 62 61 !! * Substitutions 63 # include "domzgr_substitute.h90"64 62 # include "vectopt_loop_substitute.h90" 65 63 !!---------------------------------------------------------------------- … … 118 116 DO jj = 1, jpj 119 117 DO ji = 1, jpi 120 zsfc = e1t(ji,jj) * fse3t(ji,jj,jk)118 zsfc = e1t(ji,jj) * e3t_n(ji,jj,jk) 121 119 zmask(ji,jj,jk) = tmask(ji,jj,jk) * zsfc 122 120 zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc 123 121 zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc 124 END DO125 END DO126 END DO122 END DO 123 END DO 124 END DO 127 125 DO jn = 1, nptr 128 126 zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) ) -
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r5930 r6140 30 30 USE zdf_oce ! ocean vertical physics 31 31 USE ldftra ! lateral physics: eddy diffusivity coef. 32 USE ldfdyn ! lateral physics: eddy viscosity coef. 32 33 USE sbc_oce ! Surface boundary condition: ocean fields 33 34 USE sbc_ice ! Surface boundary condition: ice fields … … 40 41 USE zdfddm ! vertical physics: double diffusion 41 42 USE diahth ! thermocline diagnostics 43 ! 42 44 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 43 45 USE in_out_manager ! I/O manager 44 USE diadimg ! dimg direct access file format output 46 USE diatmb ! Top,middle,bottom output 47 USE dia25h ! 25h Mean output 45 48 USE iom 46 49 USE ioipsl … … 53 56 USE lib_mpp ! MPP library 54 57 USE timing ! preformance summary 58 USE diurnal_bulk ! diurnal warm layer 59 USE cool_skin ! Cool skin 55 60 USE wrk_nemo ! working array 56 61 … … 74 79 !! * Substitutions 75 80 # include "zdfddm_substitute.h90" 76 # include "domzgr_substitute.h90"77 81 # include "vectopt_loop_substitute.h90" 78 82 !!---------------------------------------------------------------------- … … 97 101 END FUNCTION dia_wri_alloc 98 102 99 #if defined key_dimgout100 !!----------------------------------------------------------------------101 !! 'key_dimgout' DIMG output file102 !!----------------------------------------------------------------------103 # include "diawri_dimg.h90"104 105 #else106 103 !!---------------------------------------------------------------------- 107 104 !! Default option NetCDF output file 108 105 !!---------------------------------------------------------------------- 109 # 106 #if defined key_iomput 110 107 !!---------------------------------------------------------------------- 111 108 !! 'key_iomput' use IOM library … … 143 140 ENDIF 144 141 145 IF( .NOT.lk_vvl) THEN146 CALL iom_put( "e3t" , fse3t_n(:,:,:) )147 CALL iom_put( "e3u" , fse3u_n(:,:,:) )148 CALL iom_put( "e3v" , fse3v_n(:,:,:) )149 CALL iom_put( "e3w" , fse3w_n(:,:,:) )142 IF( ln_linssh ) THEN 143 CALL iom_put( "e3t" , e3t_n(:,:,:) ) 144 CALL iom_put( "e3u" , e3u_n(:,:,:) ) 145 CALL iom_put( "e3v" , e3v_n(:,:,:) ) 146 CALL iom_put( "e3w" , e3w_n(:,:,:) ) 150 147 ENDIF 151 148 … … 204 201 CALL iom_put( "sbu", z2d ) ! bottom i-current 205 202 ENDIF 206 207 IF ( ln_dynspg_ts ) THEN208 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current209 ELSE210 CALL iom_put( "ubar", un_b(:,:) ) ! barotropic i-current211 ENDIF212 203 213 204 CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current … … 223 214 ENDIF 224 215 225 IF ( ln_dynspg_ts ) THEN226 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic j-current227 ELSE228 CALL iom_put( "vbar", vn_b(:,:) ) ! barotropic j-current229 ENDIF230 231 216 CALL iom_put( "woce", wn ) ! vertical velocity 232 217 IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value … … 266 251 DO jj = 1, jpj 267 252 DO ji = 1, jpi 268 z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)253 z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 269 254 END DO 270 255 END DO … … 278 263 DO jj = 1, jpj 279 264 DO ji = 1, jpi 280 z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)265 z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 281 266 END DO 282 267 END DO … … 290 275 DO jj = 2, jpjm1 291 276 DO ji = fs_2, fs_jpim1 ! vector opt. 292 zztmp = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )293 zztmpx = 0.5 * ( un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) &294 & + un(ji ,jj,jk) * un(ji ,jj,jk) * e2u(ji ,jj) * fse3u(ji ,jj,jk) ) &277 zztmp = 1._wp / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 278 zztmpx = 0.5 * ( un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) & 279 & + un(ji ,jj,jk) * un(ji ,jj,jk) * e2u(ji ,jj) * e3u_n(ji ,jj,jk) ) & 295 280 & * zztmp 296 281 ! 297 zztmpy = 0.5 * ( vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) &298 & + vn(ji,jj ,jk) * vn(ji,jj ,jk) * e1v(ji,jj ) * fse3v(ji,jj ,jk) ) &282 zztmpy = 0.5 * ( vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) & 283 & + vn(ji,jj ,jk) * vn(ji,jj ,jk) * e1v(ji,jj ) * e3v_n(ji,jj ,jk) ) & 299 284 & * zztmp 300 285 ! … … 311 296 z3d(:,:,jpk) = 0.e0 312 297 DO jk = 1, jpkm1 313 z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)298 z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk) 314 299 END DO 315 300 CALL iom_put( "u_masstr", z3d ) ! mass transport in i-direction … … 346 331 z3d(:,:,jpk) = 0.e0 347 332 DO jk = 1, jpkm1 348 z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)333 z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk) 349 334 END DO 350 335 CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction … … 380 365 CALL wrk_dealloc( jpi , jpj, jpk , z3d ) 381 366 ! 367 ! If we want tmb values 368 369 IF (ln_diatmb) THEN 370 CALL dia_tmb 371 ENDIF 372 IF (ln_dia25h) THEN 373 CALL dia_25h( kt ) 374 ENDIF 375 382 376 IF( nn_timing == 1 ) CALL timing_stop('dia_wri') 383 377 ! … … 410 404 INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers 411 405 INTEGER :: jn, ierror ! local integers 412 REAL(wp) :: zsto, zout, zmax, zjulian , zdt! local scalars406 REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars 413 407 ! 414 408 REAL(wp), POINTER, DIMENSION(:,:) :: zw2d ! 2D workspace … … 418 412 IF( nn_timing == 1 ) CALL timing_start('dia_wri') 419 413 ! 420 CALL wrk_alloc( jpi,jpj , zw2d )421 IF( lk_vvl) CALL wrk_alloc( jpi,jpj,jpk , zw3d )414 CALL wrk_alloc( jpi,jpj , zw2d ) 415 IF( .NOT.ln_linssh ) CALL wrk_alloc( jpi,jpj,jpk , zw3d ) 422 416 ! 423 417 ! Output the initial state and forcings … … 435 429 436 430 ! Define frequency of output and means 437 zdt = rdt438 IF( nacc == 1 ) zdt = rdtmin439 431 clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes) 440 432 #if defined key_diainstant 441 zsto = nwrite * zdt433 zsto = nwrite * rdt 442 434 clop = "inst("//TRIM(clop)//")" 443 435 #else 444 zsto= zdt436 zsto=rdt 445 437 clop = "ave("//TRIM(clop)//")" 446 438 #endif 447 zout = nwrite * zdt448 zmax = ( nitend - nit000 + 1 ) * zdt439 zout = nwrite * rdt 440 zmax = ( nitend - nit000 + 1 ) * rdt 449 441 450 442 ! Define indices of the horizontal output zoom and vertical limit storage … … 488 480 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit 489 481 & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & 490 & nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )482 & nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set ) 491 483 CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept 492 484 & "m", ipk, gdept_1d, nz_T, "down" ) … … 524 516 CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu 525 517 & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & 526 & nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )518 & nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set ) 527 519 CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept 528 520 & "m", ipk, gdept_1d, nz_U, "down" ) … … 537 529 CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv 538 530 & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & 539 & nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )531 & nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set ) 540 532 CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept 541 533 & "m", ipk, gdept_1d, nz_V, "down" ) … … 550 542 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit 551 543 & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & 552 & nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )544 & nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set ) 553 545 CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw 554 546 & "m", ipk, gdepw_1d, nz_W, "down" ) … … 562 554 CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn 563 555 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 564 IF( lk_vvl) THEN556 IF( .NOT.ln_linssh ) THEN 565 557 CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n 566 558 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) … … 583 575 CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx 584 576 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 585 IF( .NOT. lk_vvl) THEN577 IF( ln_linssh ) THEN 586 578 CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem) 587 579 & , "KgC/m2/s", & ! sosst_cd … … 729 721 ENDIF 730 722 731 IF( lk_vvl) THEN732 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content733 CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content734 CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content735 CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content723 IF( .NOT.ln_linssh ) THEN 724 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content 725 CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content 726 CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content 727 CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content 736 728 ELSE 737 729 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature … … 740 732 CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity 741 733 ENDIF 742 IF( lk_vvl) THEN743 zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2744 CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness745 CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth734 IF( .NOT.ln_linssh ) THEN 735 zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 736 CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness 737 CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth 746 738 CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation 747 739 ENDIF … … 752 744 ! (includes virtual salt flux beneath ice 753 745 ! in linear free surface case) 754 IF( .NOT. lk_vvl) THEN746 IF( ln_linssh ) THEN 755 747 zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem) 756 748 CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst … … 837 829 ENDIF 838 830 ! 839 CALL wrk_dealloc( jpi , jpj , zw2d )840 IF( lk_vvl) CALL wrk_dealloc( jpi , jpj , jpk , zw3d )831 CALL wrk_dealloc( jpi , jpj , zw2d ) 832 IF( .NOT.ln_linssh ) CALL wrk_dealloc( jpi , jpj , jpk , zw3d ) 841 833 ! 842 834 IF( nn_timing == 1 ) CALL timing_stop('dia_wri') 843 835 ! 844 836 END SUBROUTINE dia_wri 845 # endif846 847 837 #endif 848 838 … … 867 857 INTEGER :: id_i , nz_i, nh_i 868 858 INTEGER, DIMENSION(1) :: idex ! local workspace 869 REAL(wp) :: zsto, zout, zmax, zjulian , zdt859 REAL(wp) :: zsto, zout, zmax, zjulian 870 860 !!---------------------------------------------------------------------- 871 861 ! … … 876 866 clname = cdfile_name 877 867 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname) 878 zdt = rdt879 868 zsto = rdt 880 869 clop = "inst(x)" ! no use of the mask value (require less cpu time) 881 870 zout = rdt 882 zmax = ( nitend - nit000 + 1 ) * zdt871 zmax = ( nitend - nit000 + 1 ) * rdt 883 872 884 873 IF(lwp) WRITE(numout,*) … … 895 884 zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment 896 885 CALL histbeg( clname, jpi, glamt, jpj, gphit, & 897 1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit886 1, jpi, 1, jpj, nit000-1, zjulian, rdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit 898 887 CALL histvert( id_i, "deptht", "Vertical T levels", & ! Vertical grid : gdept 899 888 "m", jpk, gdept_1d, nz_i, "down") … … 913 902 CALL histdef( id_i, "vovecrtz", "Vertical Velocity" , "m/s" , & ! vertical current 914 903 & jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 904 ! 905 CALL histdef( id_i, "ahtu" , "u-eddy diffusivity" , "m2/s" , & ! zonal current 906 & jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 907 CALL histdef( id_i, "ahtv" , "v-eddy diffusivity" , "m2/s" , & ! meridonal current 908 & jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 909 CALL histdef( id_i, "ahmt" , "t-eddy viscosity" , "m2/s" , & ! zonal current 910 & jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 911 CALL histdef( id_i, "ahmf" , "f-eddy viscosity" , "m2/s" , & ! meridonal current 912 & jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 913 ! 915 914 CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S", & ! net freshwater 916 915 & jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) … … 925 924 CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2" , & ! j-wind stress 926 925 & jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 927 IF( lk_vvl ) THEN 928 CALL histdef( id_i, "vovvldep", "T point depth" , "m" , & ! t-point depth 926 IF( .NOT.ln_linssh ) THEN 927 CALL histdef( id_i, "vovvldep", "T point depth" , "m" , & ! t-point depth 928 & jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 929 CALL histdef( id_i, "vovvle3t", "T point thickness" , "m" , & ! t-point depth 929 930 & jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 930 931 ENDIF … … 952 953 CALL histwrite( id_i, "vomecrty", kt, vn , jpi*jpj*jpk, idex ) ! now j-velocity 953 954 CALL histwrite( id_i, "vovecrtz", kt, wn , jpi*jpj*jpk, idex ) ! now k-velocity 955 ! 956 CALL histwrite( id_i, "ahtu" , kt, ahtu , jpi*jpj*jpk, idex ) ! aht at u-point 957 CALL histwrite( id_i, "ahtv" , kt, ahtv , jpi*jpj*jpk, idex ) ! - at v-point 958 CALL histwrite( id_i, "ahmt" , kt, ahmt , jpi*jpj*jpk, idex ) ! ahm at t-point 959 CALL histwrite( id_i, "ahmf" , kt, ahmf , jpi*jpj*jpk, idex ) ! - at f-point 960 ! 954 961 CALL histwrite( id_i, "sowaflup", kt, emp-rnf , jpi*jpj , idex ) ! freshwater budget 955 962 CALL histwrite( id_i, "sohefldo", kt, qsr + qns , jpi*jpj , idex ) ! total heat flux … … 959 966 CALL histwrite( id_i, "sometauy", kt, vtau , jpi*jpj , idex ) ! j-wind stress 960 967 968 IF( .NOT.ln_linssh ) THEN 969 CALL histwrite( id_i, "vovvldep", kt, gdept_n(:,:,:), jpi*jpj*jpk, idex )! T-cell depth 970 CALL histwrite( id_i, "vovvle3t", kt, e3t_n (:,:,:) , jpi*jpj*jpk, idex )! T-cell thickness 971 END IF 961 972 ! 3. Close the file 962 973 ! ----------------- 963 974 CALL histclo( id_i ) 964 #if ! defined key_iomput && ! defined key_dimgout975 #if ! defined key_iomput 965 976 IF( ninist /= 1 ) THEN 966 977 CALL histclo( nid_T ) … … 972 983 ! 973 984 END SUBROUTINE dia_wri_state 985 974 986 !!====================================================================== 975 987 END MODULE diawri
Note: See TracChangeset
for help on using the changeset viewer.
