source: trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90 @ 848

Last change on this file since 848 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

  • nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this
  • Property svn:executable set to *
File size: 10.1 KB
Line 
1subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf,          &
2     qsurf,dqsurf,dqs_hyd,pcapcal,                &
3     albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)
4
5  use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
6  USE surfdat_h
7  use comdiurn_h
8  USE comgeomfi_h
9  USE tracer_h
10
11  implicit none
12
13!==================================================================
14!     
15!     Purpose
16!     -------
17!     Calculate the surface hydrology and albedo changes.
18!     
19!     Authors
20!     -------
21!     Adapted from LMDTERRE by B. Charnay (2010). Further
22!     modifications by R. Wordsworth (2010).
23!     
24!     Called by
25!     ---------
26!     physiq.F
27!     
28!     Calls
29!     -----
30!     none
31!
32!     Notes
33!     -----
34!     rnat is terrain type: 0-ocean; 1-continent
35!     
36!==================================================================
37
38#include "dimensions.h"
39#include "dimphys.h"
40#include "comcstfi.h"
41#include "callkeys.h"
42
43      integer ngrid,nq
44
45!     Inputs
46!     ------
47      real albedoocean
48      parameter (albedoocean=0.07)
49      real albedoice
50      save albedoice
51
52      real snowlayer
53      parameter (snowlayer=33.0)        ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
54      real oceantime
55      parameter (oceantime=10*24*3600)
56
57      logical oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
58      logical activerunoff ! enable simple runoff scheme?
59      logical oceanalbvary ! ocean albedo varies with the diurnal cycle?
60
61!     Arguments
62!     ---------
63      integer rnat(ngrid) ! I changed this to integer (RW)
64      real,dimension(:),allocatable,save :: runoff
65      real totalrunoff, tsea, oceanarea
66      save oceanarea
67
68      real ptimestep
69      real mu0(ngrid)
70      real qsurf(ngrid,nq), tsurf(ngrid)
71      real dqsurf(ngrid,nq), pdtsurf(ngrid)
72      real hice(ngrid)
73      real albedo0(ngrid), albedo(ngrid)
74
75      real oceanarea2
76
77!     Output
78!     ------
79      real dqs_hyd(ngrid,nq)
80      real pdtsurf_hyd(ngrid)
81
82!     Local
83!     -----
84      real a,b,E
85      integer ig,iq, icap ! wld like to remove icap
86      real fsnoi, subli, fauxo
87      real twater(ngrid)
88      real pcapcal(ngrid)
89      real hicebis(ngrid)
90      real zqsurf(ngrid,nq)
91      real ztsurf(ngrid)
92
93      integer ivap, iliq, iice
94      save ivap, iliq, iice
95
96      logical firstcall
97      save firstcall
98
99      data firstcall /.true./
100
101      oceanbulkavg=.false.
102      activerunoff=.false.
103      oceanalbvary=.false.
104
105      if(firstcall)then
106
107         ALLOCATE(runoff(ngrid))
108
109         ivap=igcm_h2o_vap
110         iliq=igcm_h2o_vap
111         iice=igcm_h2o_ice
112       
113         write(*,*) "hydrol: ivap=",ivap
114         write(*,*) "        iliq=",iliq
115         write(*,*) "        iice=",iice
116
117!     Here's the deal: iice is used in place of igcm_h2o_ice both on the
118!                      surface and in the atmosphere. ivap is used in
119!                      place of igcm_h2o_vap ONLY in the atmosphere, while
120!                      iliq is used in place of igcm_h2o_vap ONLY on the
121!                      surface.
122
123!                      Soon to be extended to the entire water cycle...
124
125!     Ice albedo = snow albedo for now
126         albedoice=albedosnow
127
128!     Total ocean surface area
129         oceanarea=0.
130         do ig=1,ngrid
131            if(rnat(ig).eq.0)then
132               oceanarea=oceanarea+area(ig)
133            endif
134         enddo
135
136         if(oceanbulkavg.and.(oceanarea.le.0.))then
137            print*,'How are we supposed to average the ocean'
138            print*,'temperature, when there are no oceans?'
139            call abort
140         endif
141
142         if(activerunoff.and.(oceanarea.le.0.))then
143            print*,'You have enabled runoff, but you have no oceans.'
144            print*,'Where did you think the water was going to go?'
145            call abort
146         endif
147         
148         firstcall = .false.
149      endif
150
151!     add physical tendencies already calculated
152!     ------------------------------------------
153
154      do ig=1,ngrid
155         ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
156         pdtsurf_hyd(ig)=0.0
157         do iq=1,nq
158            zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
159         enddo   
160      enddo
161 
162      do ig=1,ngrid
163         do iq=1,nq
164            dqs_hyd(ig,iq) = 0.0
165         enddo
166      enddo
167
168      do ig = 1, ngrid
169
170!     Ocean
171!     -----
172         if(rnat(ig).eq.0)then
173
174!     re-calculate oceanic albedo
175            if(diurnal.and.oceanalbvary)then
176               fauxo      = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
177               albedo(ig) = 1.1*( .03 + .630/( 1. + fauxo*fauxo))
178               albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
179            else
180               albedo(ig) = albedoocean ! modif Benjamin
181            end if
182
183!     calculate oceanic ice height including the latent heat of ice formation
184!     hice is the height of oceanic ice with a maximum of maxicethick.
185            hice(ig)    = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
186
187!            twater(ig)  = tsurf(ig) + ptimestep*zdtsurf(ig) &
188            twater(ig)  = ztsurf(ig) - hice(ig)*RLFTT*rhowater/pcapcal(ig)
189            ! this is temperature water would have if we melted entire ocean ice layer
190            hicebis(ig) = hice(ig)
191            hice(ig)    = 0.
192
193            if(twater(ig) .lt. T_h2O_ice_liq)then
194               E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick)
195               hice(ig)        = E/(RLFTT*rhowater)
196               hice(ig)        = max(hice(ig),0.0)
197               hice(ig)        = min(hice(ig),maxicethick)
198               pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT*rhowater/pcapcal(ig)/ptimestep             
199               albedo(ig)      = albedoice
200
201!               if (zqsurf(ig,iice).ge.snowlayer) then
202!                  albedo(ig) = albedoice
203!               else
204!                  albedo(ig) = albedoocean &
205!                       + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
206!               endif
207
208            else
209
210               pdtsurf_hyd(ig) = -hicebis(ig)*RLFTT*rhowater/pcapcal(ig)/ptimestep               
211               albedo(ig)      = albedoocean
212
213            endif
214
215            zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
216            zqsurf(ig,iice) = hice(ig)*rhowater
217
218
219!     Continent
220!     ---------
221         elseif (rnat(ig).eq.1) then
222
223!     melt the snow
224            if(ztsurf(ig).gt.T_h2O_ice_liq)then
225               if(zqsurf(ig,iice).gt.1.0e-8)then
226
227                  a     = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
228                  b     = zqsurf(ig,iice)
229                  fsnoi = min(a,b)
230
231                  zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi
232                  zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi
233
234!                 thermal effects
235                  pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep 
236
237               endif
238            else
239
240!     freeze the water
241               if(zqsurf(ig,iliq).gt.1.0e-8)then
242
243                  a     = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
244                  b     = zqsurf(ig,iliq)
245                 
246                  fsnoi = min(a,b)
247
248                  zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi
249                  zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi
250
251!                 thermal effects
252                  pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep 
253
254               endif
255            endif
256           
257!     deal with runoff
258            if(activerunoff)then
259
260               runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
261               if(ngrid.gt.1)then ! runoff only exists in 3D
262                  if(runoff(ig).ne.0.0)then
263                     zqsurf(ig,iliq) = mx_eau_sol
264!                    runoff is added to ocean at end
265                  endif
266               end if
267
268            endif
269
270!     re-calculate continental albedo
271            albedo(ig) = albedo0(ig) ! albedo0 = base values
272            if (zqsurf(ig,iice).ge.snowlayer) then
273               albedo(ig) = albedosnow
274            else
275               albedo(ig) = albedo0(ig) &
276                    + (albedosnow - albedo0(ig))*zqsurf(ig,iice)/snowlayer
277            endif
278
279         else
280
281            print*,'Surface type not recognised in hydrol.F!'
282            print*,'Exiting...'
283            call abort
284
285         endif
286
287      end do ! ig=1,ngrid
288
289!     perform crude bulk averaging of temperature in ocean
290!     ----------------------------------------------------
291      if(oceanbulkavg)then
292
293         oceanarea2=0.       
294         DO ig=1,ngrid
295            if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
296               oceanarea2=oceanarea2+area(ig)*pcapcal(ig)
297            end if
298         END DO
299       
300         tsea=0.
301         DO ig=1,ngrid
302            if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then       
303               tsea=tsea+ztsurf(ig)*area(ig)*pcapcal(ig)/oceanarea2
304            end if
305         END DO
306
307         DO ig=1,ngrid
308            if((rnat(ig).eq.0).and.(hice(ig).eq.0))then
309               pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
310            end if
311         END DO
312
313         print*,'Mean ocean temperature               = ',tsea,' K'
314
315      endif
316
317!     shove all the runoff water into the ocean
318!     -----------------------------------------
319      if(activerunoff)then
320
321         totalrunoff=0.
322         do ig=1,ngrid
323            if (rnat(ig).eq.1) then
324               totalrunoff = totalrunoff + area(ig)*runoff(ig)
325            endif
326         enddo
327         
328         do ig=1,ngrid
329            if (rnat(ig).eq.0) then
330               zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
331                    totalrunoff/oceanarea
332            endif
333         enddo
334
335      endif
336
337!     Re-add the albedo effects of CO2 ice if necessary
338!     -------------------------------------------------
339      if(co2cond)then
340
341         icap=1
342         do ig=1,ngrid
343            if (qsurf(ig,igcm_co2_ice).gt.0) then
344               albedo(ig) = albedice(icap)
345            endif
346         enddo
347         
348      endif
349
350
351      do ig=1,ngrid
352         dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep
353         dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep
354      enddo
355
356      call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
357
358      return
359    end subroutine hydrol
Note: See TracBrowser for help on using the repository browser.