source: LMDZ5/trunk/libf/phylmd/surf_ocean_mod.F90 @ 2390

Last change on this file since 2390 was 2322, checked in by lguez, 9 years ago

Bug fix in alboc_cd (more than 11 years old). rmu0 must not be
modified in alboc_cd. The corresponding actual argument in subroutine
surf_ocean is an intent(in) argument of surf_ocean. It probably worked
before because the compiler made a copy of the argument but I would
say the behavior was in principle compiler-dependent and
optimization-level-dependent.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 9.1 KB
Line 
1!
2MODULE surf_ocean_mod
3
4  IMPLICIT NONE
5
6CONTAINS
7!
8!******************************************************************************
9!
10  SUBROUTINE surf_ocean(rlon, rlat, swnet, lwnet, alb1, &
11       windsp, rmu0, fder, tsurf_in, &
12       itime, dtime, jour, knon, knindex, &
13       p1lay, z1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
14       AcoefH, AcoefQ, BcoefH, BcoefQ, &
15       AcoefU, AcoefV, BcoefU, BcoefV, &
16       ps, u1, v1, gustiness, rugoro, pctsrf, &
17       snow, qsurf, agesno, &
18       z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, &
19       tsurf_new, dflux_s, dflux_l, lmt_bils, &
20       flux_u1, flux_v1)
21
22  use albedo, only: alboc, alboc_cd
23  USE dimphy
24  USE surface_data, ONLY     : type_ocean
25  USE ocean_forced_mod, ONLY : ocean_forced_noice
26  USE ocean_slab_mod, ONLY   : ocean_slab_noice
27  USE ocean_cpl_mod, ONLY    : ocean_cpl_noice
28  USE indice_sol_mod
29!
30! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force,
31! slab or couple). The calculations of albedo and rugosity for the ocean surface are
32! done in here because they are identical for the different modes of ocean.
33
34
35    INCLUDE "YOMCST.h"
36
37    include "clesphys.h"
38    ! for cycle_diurne
39
40! Input variables
41!******************************************************************************
42    INTEGER, INTENT(IN)                      :: itime, jour, knon
43    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
44    REAL, INTENT(IN)                         :: dtime
45    REAL, DIMENSION(klon), INTENT(IN)        :: rlon, rlat
46    REAL, DIMENSION(klon), INTENT(IN)        :: swnet  ! net shortwave radiation at surface 
47    REAL, DIMENSION(klon), INTENT(IN)        :: lwnet  ! net longwave radiation at surface 
48    REAL, DIMENSION(klon), INTENT(IN)        :: alb1   ! albedo in visible SW interval
49    REAL, DIMENSION(klon), INTENT(IN)        :: windsp
50    REAL, DIMENSION(klon), INTENT(IN)        :: rmu0 
51    REAL, DIMENSION(klon), INTENT(IN)        :: fder
52    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
53    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay,z1lay ! pression (Pa) et altitude (m) du premier niveau
54    REAL, DIMENSION(klon), INTENT(IN)        :: cdragh
55    REAL, DIMENSION(klon), INTENT(IN)        :: cdragm
56    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
57    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
58    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefH, AcoefQ, BcoefH, BcoefQ
59    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefU, AcoefV, BcoefU, BcoefV
60    REAL, DIMENSION(klon), INTENT(IN)        :: ps
61    REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
62    REAL, DIMENSION(klon), INTENT(IN)        :: rugoro
63    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
64
65! In/Output variables
66!******************************************************************************
67    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
68    REAL, DIMENSION(klon), INTENT(INOUT)     :: qsurf
69    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno
70
71! Output variables
72!******************************************************************************
73    REAL, DIMENSION(klon), INTENT(OUT)       :: z0m, z0h
74!albedo SB >>>
75!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new  ! new albedo in visible SW interval
76!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new  ! new albedo in near IR interval
77    REAL, DIMENSION(6), INTENT(IN)          :: SFRWL
78    REAL, DIMENSION(klon,nsw), INTENT(OUT)       :: alb_dir_new,alb_dif_new
79!albedo SB <<<     
80    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
81    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
82    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l     
83    REAL, DIMENSION(klon), INTENT(OUT)       :: lmt_bils
84    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
85
86! Local variables
87!******************************************************************************
88    INTEGER               :: i, k
89    REAL                  :: tmp
90    REAL, PARAMETER       :: cepdu2=(0.1)**2
91    REAL, DIMENSION(klon) :: alb_eau
92    REAL, DIMENSION(klon) :: radsol
93    REAL, DIMENSION(klon) :: cdragq ! Cdrag pour l'evaporation
94
95! End definition
96!******************************************************************************
97
98
99!******************************************************************************
100! Calculate total net radiance at surface
101!
102!******************************************************************************
103    radsol(:) = 0.0
104    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
105
106!******************************************************************************
107! Cdragq computed from cdrag
108! The difference comes only from a factor (f_z0qh_oce) on z0, so that
109! it can be computed inside surf_ocean
110! More complicated appraches may require the propagation through
111! pbl_surface of an independant cdragq variable.
112!******************************************************************************
113
114    IF ( f_z0qh_oce .ne. 1.) THEN
115! Si on suit les formulations par exemple de Tessel, on
116! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55
117       cdragq(:)=cdragh(:)*                                      &
118       log(z1lay(:)/z0h(:))/log(z1lay(:)/(f_z0qh_oce*z0h(:)))
119    ELSE
120       cdragq(:)=cdragh(:)
121    ENDIF
122
123!******************************************************************************
124! Switch according to type of ocean (couple, slab or forced)
125!******************************************************************************
126    SELECT CASE(type_ocean)
127    CASE('couple')
128       CALL ocean_cpl_noice( &
129            swnet, lwnet, alb1, &
130            windsp, fder, &
131            itime, dtime, knon, knindex, &
132            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow,temp_air,spechum,&
133            AcoefH, AcoefQ, BcoefH, BcoefQ, &
134            AcoefU, AcoefV, BcoefU, BcoefV, &
135            ps, u1, v1, gustiness, &
136            radsol, snow, agesno, &
137            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
138            tsurf_new, dflux_s, dflux_l)
139
140    CASE('slab')
141       CALL ocean_slab_noice( &
142            itime, dtime, jour, knon, knindex, &
143            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum,&
144            AcoefH, AcoefQ, BcoefH, BcoefQ, &
145            AcoefU, AcoefV, BcoefU, BcoefV, &
146            ps, u1, v1, gustiness, tsurf_in, &
147            radsol, snow, &
148            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
149            tsurf_new, dflux_s, dflux_l, lmt_bils)
150       
151    CASE('force')
152       CALL ocean_forced_noice( &
153            itime, dtime, jour, knon, knindex, &
154            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, &
155            temp_air, spechum, &
156            AcoefH, AcoefQ, BcoefH, BcoefQ, &
157            AcoefU, AcoefV, BcoefU, BcoefV, &
158            ps, u1, v1, gustiness, &
159            radsol, snow, agesno, &
160            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
161            tsurf_new, dflux_s, dflux_l)
162    END SELECT
163
164!******************************************************************************
165! fcodron: compute lmt_bils  forced case (same as wfbils_oce / 1.-contfracatm)
166!******************************************************************************
167    IF (type_ocean.NE.'slab') THEN
168        lmt_bils(:)=0.
169        DO i=1,knon
170           lmt_bils(knindex(i))=(swnet(i)+lwnet(i)+fluxsens(i)+fluxlat(i)) &
171           *pctsrf(knindex(i),is_oce)/(1.-zmasq(knindex(i)))
172        END DO
173    END IF
174
175!******************************************************************************
176! Calculate albedo
177!******************************************************************************
178!albedo SB >>>
179  if(iflag_albedo==1)then
180    call ocean_albedo(knon,rmu0,knindex,windsp,SFRWL,alb_dir_new,alb_dif_new)
181  else
182    IF (cycle_diurne) THEN
183       CALL alboc_cd(rmu0,alb_eau)
184    ELSE
185       CALL alboc(REAL(jour),rlat,alb_eau)
186    ENDIF
187
188    DO i =1, knon
189      do  k=1,nsw
190       alb_dir_new(i,k) = alb_eau(knindex(i))
191      enddo
192    ENDDO
193     alb_dif_new=0.05 !alb_dir_new
194endif
195
196!albedo SB <<<
197
198!******************************************************************************
199! Calculate the rugosity
200!******************************************************************************
201IF (iflag_z0_oce==0) THEN
202    DO i = 1, knon
203       tmp = MAX(cepdu2,gustiness(i)+u1(i)**2+v1(i)**2)
204       z0m(i) = 0.018*cdragm(i) * (gustiness(i)+u1(i)**2+v1(i)**2)/RG  &
205            +  0.11*14e-6 / SQRT(cdragm(i) * tmp)
206       z0m(i) = MAX(1.5e-05,z0m(i))
207    ENDDO   
208    z0h(1:knon)=z0m(1:knon) ! En attendant mieux
209
210ELSE IF (iflag_z0_oce==1) THEN
211    DO i = 1, knon
212       tmp = MAX(cepdu2,gustiness(i)+u1(i)**2+v1(i)**2)
213       z0m(i) = 0.018*cdragm(i) * (gustiness(i)+u1(i)**2+v1(i)**2)/RG  &
214            + 0.11*14e-6 / SQRT(cdragm(i) * tmp)
215       z0m(i) = MAX(1.5e-05,z0m(i))
216       z0h(i)=0.4*14e-6 / SQRT(cdragm(i) * tmp)
217    ENDDO
218ELSE
219       STOP'version non prevue'
220ENDIF
221!
222!******************************************************************************
223  END SUBROUTINE surf_ocean
224!******************************************************************************
225!
226END MODULE surf_ocean_mod
Note: See TracBrowser for help on using the repository browser.